Building An Restrained Nonbonded Model

A restrained nonbonded model could be used when the normal bonded model failed under certain circumistances. For example, for the case that there are two water molecules connected to the same metal ion, which would cause structure distortion because of the 1-4 VDW interaction between oxygen in one water and hydrogen in the other is zero while the 1-4 electrostatic interaction between them are nonzero and attractive. The strategy is to delete related bonds in the LEaP modeling, treat the bonds in a nonbonded way, while use harmonic restraints for these bonds in the minimization and molecular dynamics simulations. In this way the full nonbonded interactions between two water molecules would be considered to solve the problem. Users can treat all the coordination bonds in this way or treat some of them in the bonded way while the others in this way (a hybrid bonded/restrained nonbonded model). Here we show one example based on the example illustrated in the MCPB.py tutorial:

1. Detect the bond information related to the metal site:

Here we can still use the "printBonds :ZN1" command in ParmEd to get the metal ion related bond information. We can get the following information:

2. Generate new topology and coordinate files (effort: several minutes):

Here we delete certain bond command(s) in LEaP input file and rename the output pdb, prmtop and inpcrd files. Here we delete the bonds between zinc ion and all the ligating atoms, treat all the four coordinate bonds with the nonbonded model. Users may use a hybrid model with only delete specific bond command(s) in the tleap file based on their own choice. Reminder one only needs to apply harmonic restraint(s) to the bond(s) deleted in the LEaP modeling during the minimization and molecular dynamics simualtion. Here is the tleap input file we use: (1OKL_rnb_tleap.in). We perform the command:

tleap -s -f 1OKL_rnb_tleap.in > 1OKL_rnb_tleap.out

Now we have the new topology and coordinate files 1OKL_rnb_solv.prmtop and 1OKL_rnb_solv.inpcrd generated.

3. Apply the harmonic restraints to the bonds:

Here we created an RST file for the harmonic restraints: 1OKL_rnb.RST remind that the strength of the harmonic restraints should be consistent with the Figure shown above:

# Harmonic restraints for the restrained nonbonded model of 1OKL &rst iat=1401,4018, r1=0., r2=2.0229, r3=2.0229, r4=100., rk2=72.7, rk3=72.7,/ &rst iat=1438,4018, r1=0., r2=2.0610, r3=2.0610, r4=100., rk2=57.6, rk3=57.6,/ &rst iat=1782,4018, r1=0., r2=2.0354, r3=2.0354, r4=100., rk2=66.3, rk3=66.3,/ &rst iat=4018,4035, r1=0., r2=1.9619, r3=1.9619, r4=100., rk2=101.2, rk3=101.2,/

Herein iat means the atomic IDs of two atoms in the coordinate bond; r2 and r3 correspond to the bond length of the coordinate bond; r1 and r4 are the lower and upper limits of the bond length to apply the restraints; rk2 and rk3 correspond to the force constant of the coordinate bond. This kind of restraints was also used in the umbrella sampling, one can check the example in Page 420 of the AMBER 2017 manual for how to apply these restraints into the minimization or simulation. One can also go to the NMR refinement tutorial to be more familiar with these restraints.