AMBER Archive (2006)

Subject: RE: AMBER: how to restrain some residues in a NPT MD?

From: Ross Walker (ross_at_rosswalker.co.uk)
Date: Mon Oct 30 2006 - 10:34:46 CST

Dear Zhihong

>From your description of runs 4 and 5 below:

> 4. mpirun -np 8 pmemd.MPI -O -i equ_product.in -o
> equ_product.out -p system.prmtop -c system_heat.rst -r
> system_equ2_product.rst -x system_equ2_product.mdcrd -ref
> system_heat.rst (do NPT MD at 300K for 1ns)

and

> mpirun -np 8 pmemd.MPI -O -i equ_product.in -o product_2.out
> -p system.prmtop -c system_heat.rst -r system_product_2.rst
> -x system_product_2.mdcrd -ref system_heat.rst

I can't see what is different in the two runs - other than some output
filenames. Hence the second run should just give you the same output as the
first. However, I will assume what you actually did is:

4. mpirun -np 8 pmemd.MPI -O -i equ_product.in -o equ_product.out -p
system.prmtop -c system_heat.rst -r system_equ2_product.rst -x
system_equ2_product.mdcrd -ref system_heat.rst (do NPT MD at 300K for 1ns)

5. mpirun -np 8 pmemd.MPI -O -i equ_product.in -o product_2.out -p
system.prmtop -c system_equ2_product.rst -r system_product_2.rst -x
system_product_2.mdcrd -ref system_heat.rst

That is you used the rst file from step 4 as the input coordinates for step
5 so that you could continue the 1ns simulation for a further 1ns. If this
is not the case then I think you will have to more carefully explain what
you did. Anyway, the remainder of my answer assumes the above:

The problem you are seeing is a limitation of Amber based on the way in
which restraints work when using constant pressure. The box size can change
during a simulation and hence if you restart it with restraints pointing to
a reference structure with a different box size you can get huge restraint
forces and your system blows up.

Simply change the reference structure you use for the restarted NPT
simulation from:

mpirun -np 8 pmemd.MPI -O -i equ_product.in -o product_2.out
-p system.prmtop -c system_equ2_product.rst -r system_product_2.rst
-x system_product_2.mdcrd -ref system_heat.rst

to

mpirun -np 8 pmemd.MPI -O -i equ_product.in -o product_2.out
-p system.prmtop -c system_equ2_product.rst -r system_product_2.rst
-x system_product_2.mdcrd -ref system_equ2_product.rst

That is, the reference structure for the restraints should be the same as
the input coordinates specified by -c.

And it should be fine.

All the best
Ross

/\
\/
|\oss Walker

| HPC Consultant and Staff Scientist |
| San Diego Supercomputer Center |
| Tel: +1 858 822 0854 | EMail:- ross_at_rosswalker.co.uk |
| http://www.rosswalker.co.uk | PGP Key available on request |

Note: Electronic Mail is not secure, has no guarantee of delivery, may not
be read every day, and should not be used for urgent or sensitive issues.

> -----Original Message-----
> From: owner-amber_at_scripps.edu
> [mailto:owner-amber_at_scripps.edu] On Behalf Of Zhihong Yu
> Sent: Monday, October 30, 2006 07:21
> To: amber
> Subject: AMBER: how to restrain some residues in a NPT MD?
>
> Dear All,
>
> I have a system containing several protein fragments,
> because I extract all residues within the distance of 20A of
> an inhibitor, I want to do NPT MD with very weakly restraint
> (1.0 kcal/mol*A^2) to those residues in the 15~20A shell, so
> I did the MD in following four steps, and everything is OK.
>
> 1. mpirun -np 8 pmemd.MPI -O -i min1.in -o min1.out -p
> system.prmtop -c system.inpcrd -r system_min1.rst -ref
> system.inpcrd (restrain solute minimize water)
> 2. mpirun -np 8 pmemd.MPI -O -i min2.in -o min2.out -p
> system.prmtop -c system_min1.rst -r system_min2.rst (minimize
> the whole system)
> 3. mpirun -np 8 pmemd.MPI -O -i heat.in -o heat.out -p
> system.prmtop -c system_min2.rst -r system_heat.rst -x
> system_heat.mdcrd -ref system_min2.rst (weakly restrain
> solute and heat the system from 0K to 300K in 20ps)
> 4. mpirun -np 8 pmemd.MPI -O -i equ_product.in -o
> equ_product.out -p system.prmtop -c system_heat.rst -r
> system_equ2_product.rst -x system_equ2_product.mdcrd -ref
> system_heat.rst (do NPT MD at 300K for 1ns)
>
> When I analysis the MD trajectory, I found the backbone RMSD
> is not very stable in 1ns simulation, so I want to do another
> 1ns MD simulation based on the step 4 as follows:
>
> mpirun -np 8 pmemd.MPI -O -i equ_product.in -o product_2.out
> -p system.prmtop -c system_heat.rst -r system_product_2.rst
> -x system_product_2.mdcrd -ref system_heat.rst
>
> at this moment, I found a problem that the temperature
> increased to very high, while in step 4 the temperature is
> very stable, average temperature is 300.03K, temperature RMS
> fluctuations is 1.84K. At the same time, "Density" decreased
> to the value almost equal to the initial value in step 4(0.9571)
>
> NSTEP = 50 TIME(PS) = 1020.100 TEMP(K) =
> 483.64 PRESS = 3604.9
> -------------------
> Etot = -45928.6887 EKtot = 29393.8883 EPtot
> = -75322.5770
> BOND = 1159.7388 ANGLE = 3198.1474 DIHED
> = 3618.2022
> 1-4 NB = 1283.6419 1-4 EEL = 11045.3955 VDWAALS
> = 9122.6836
> EELEC = -111006.6701 EHBOND = 0.0000 RESTRAINT
> = 6256.2836
> EAMBER (non-restraint) = -81578.8606
> EKCMT = 15282.0460 VIRIAL = -8022.8348 VOLUME
> = 299419.6738
> Density
> = 0.9879
>
> ---------------------------
>
>
> NSTEP = 100 TIME(PS) = 1020.200 TEMP(K) =
> 483.56 PRESS = 3562.6
> -------------------
> Etot = -48449.6942 EKtot = 29388.6835 EPtot
> = -77838.3777
> BOND = 1155.9994 ANGLE = 3288.6053 DIHED
> = 3788.3022
> 1-4 NB = 1279.2409 1-4 EEL = 11012.7785 VDWAALS
> = 9584.3672
> EELEC = -109772.8083 EHBOND = 0.0000 RESTRAINT
> = 1825.1371
> EAMBER (non-restraint) = -79663.5149
> EKCMT = 15907.0210 VIRIAL = -7323.8078 VOLUME
> = 302009.3056
> Density
> = 0.9900
>
> ---------------------------
>
>
> but when I didn't restrain those residues within 15~20A shell
> and restart MD based on step 4 as follows:
>
> mpirun -np 8 pmemd.MPI -O -i equ_product.in -o product_2.out
> -p system.prmtop -c system_heat.rst -r system_product_2.rst
> -x system_product_2.mdcrd
>
> everything looks good:
>
> NSTEP = 50 TIME(PS) = 1020.100 TEMP(K) =
> 300.26 PRESS = -45.9
> Etot = -65782.2433 EKtot = 18248.4030 EPtot
> = -84030.6463
> BOND = 988.3259 ANGLE = 2784.8558 DIHED
> = 3480.3918
> 1-4 NB = 1207.0591 1-4 EEL = 11079.5310 VDWAALS
> = 7853.6265
> EELEC = -111424.4364 EHBOND = 0.0000 RESTRAINT
> = 0.0000
> EKCMT = 7287.8821 VIRIAL = 7583.4197 VOLUME
> = 298489.4102
> Density
> = 1.0137
> Ewald error estimate: 0.1396E-03
>
>
> NSTEP = 100 TIME(PS) = 1020.200 TEMP(K) =
> 301.29 PRESS = -290.3
> Etot = -65747.3646 EKtot = 18311.4818 EPtot
> = -84058.8464
> BOND = 988.4036 ANGLE = 2795.8995 DIHED
> = 3576.7810
> 1-4 NB = 1208.2661 1-4 EEL = 11100.0828 VDWAALS
> = 7787.7048
> EELEC = -111515.9842 EHBOND = 0.0000 RESTRAINT
> = 0.0000
> EKCMT = 7325.2614 VIRIAL = 9195.5766 VOLUME
> = 298360.9939
> Density
> = 1.0141
> Ewald error estimate: 0.1364E-04
>
> So I think these two problems (temperature and density) are
> related to the restrain in NPT MD, which rst file should I
> use as reference coordinates in restart NPT MD? or I can't
> make use of the finished 1ns MD and must do a new MD
> completely for 2ns as step4? Moreover it seems that it's not
> correct using "system_heat.rst" as reference coordinates in
> step4? then how to restrain some residues in a NPT MD?
>
>
> Any advice will be greatly appreciated! Thanks in advance
> & Best regards!
>
> yours sincerely, Zhihong Yu
>
> PS: "equ_product.in" file is attached for your information

-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber_at_scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo_at_scripps.edu