# 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 [1]: from prody import *
In [2]: from matplotlib.pylab import *
In [3]: ion() # turn interactive mode on
```

We start with parsing a PDB file by passing an identifier.

```
In [4]: p38 = parsePDB('1p38')
In [5]: p38
Out[5]: <AtomGroup: 1p38 (2962 atoms)>
```

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

```
In [6]: calphas = p38.select('protein and name CA')
In [7]: calphas
Out[7]: <Selection: 'protein and name CA' from 1p38 (351 atoms)>
```

## Force Constant Function¶

We define the aformentioned function as follows:

```
In [8]: 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 [9]: dist = arange(0, 20, 0.1)
In [10]: gamma = map(gammaDistanceDependent, dist ** 2)
In [11]: plot(dist, gamma, lw=4);
In [12]: axis([0, 20, 0, 12]);
In [13]: xlabel('Distance (A)');
In [14]: ylabel('Force constant');
In [15]: grid();
```

## ANM calculations¶

We use selected atoms (351 Cα’s) and `gammaDistanceDependent`

function
for ANM calculations as follows:

```
In [16]: anm = ANM('1p38')
In [17]: anm.buildHessian(calphas, cutoff=15, gamma=gammaDistanceDependent)
In [18]: anm.calcModes()
```

For more detailed examples see Anisotropic Network Model (ANM) or Gaussian Network Model (GNM).