# Custom Gamma Functions¶

This example shows how to develop custom force constant functions for `ANM` (or `GNM`) calculations.

We will use the relation shown in the figure below. For Cα atoms that are 10 to 15 Å apart from each other, we use a unit force constant. For those that are 4 to 10 Å apart, we use a 2 times stronger force constant. For those that are within 4 Å of each other (i.e. those from connected residue pairs), we use a 10 times stronger force constant.

We will obtain an `ANM` instance that stores Hessian and Kirchhoff matrices and normal mode data describing the intrinsic dynamics of the protein structure. `ANM` instances and individual normal modes (`Mode`) can be used as input to functions in `prody.dynamics` module.

## Parse structure¶

We start by importing everything from ProDy, NumPy, and Matplotlib packages:

```In : from prody import *

In : from matplotlib.pylab import *

In : ion() # turn interactive mode on
```

```In : p38 = parsePDB('1p38')

In : p38
Out: <AtomGroup: 1p38 (2962 atoms)>
```

We want to use only Cα atoms, so we select them:

```In : calphas = p38.select('protein and name CA')

In : calphas
Out: <Selection: 'protein and name CA' from 1p38 (351 atoms)>
```

## Force Constant Function¶

We define the aformentioned function as follows:

```In : def gammaDistanceDependent(dist2, *args):
...:     """Return a force constant based on the given square distance."""
...:     if dist2 <= 16:
...:         return 10
...:     elif dist2 <= 100:
...:         return 2
...:     elif dist2 <= 225:
...:         return 1
...:     else:
...:         return 0
...:
```

Note that the input to this function from `ANM` or `GNM` is the square of the distance. In addition, node (atom or residue) indices are passed to this function, that’s why we used `*args` in the function definition.

Let’s test how it works:

```In : dist = arange(0, 20, 0.1)

In : gamma = map(gammaDistanceDependent, dist ** 2)

In : plot(dist, gamma, lw=4);

In : axis([0, 20, 0, 12]);

In : xlabel('Distance (A)');

In : ylabel('Force constant');

In : grid();
``` ## ANM calculations¶

We use selected atoms (351 Cα’s) and `gammaDistanceDependent` function for ANM calculations as follows:

```In : anm = ANM('1p38')