AMBER Archive (2005)

Subject: AMBER: broblem of replica exchange MD with polarizable force field

From: wenfei Li (wfli_at_biophy.nju.edu.cn)
Date: Wed Aug 31 2005 - 19:51:32 CDT


Dear Amber users,

I'm doing a replica-exchange MD simulation with the polarizable force field (ff02). In calculating the induced dipole, the Car-Parinello scheme is used. However, the polarization part of the potential energy is always positive which is obviously unreasonable. Meanwhile, the polarization part of the potential energy is zero at the first time step of the MD. After reading the code carefully, I found that the varialble "initdip" which appears at Line 64 of ew_handle_dips.f is always zero after the first attempt of exchange:

    "if ( indmeth < 3 .or. initdip == 1 )then)."

This results in that the initial value of the induced dipoles are always zero. I guess that the simultaneously use of the multisander and "data initdip/1/ " at Line 39 of ew_handle_dips.f may result in the above problem.

I tried to change the Line 64 of ew_handle_dips.f from "if ( indmeth < 3 .or. initdip == 1 )then)." to "if ( indmeth < 3 .or. nstep == 0 )then)."and added a common block "common /nstep_dip/ nstep" to both runmd.f and ew_handle_dips.f. After the above modifications,the results seems reasonable(The polarization part of the potential energy is negative). Is this a bug or something wrong with my input data?

Here is one of my input data:

Zinc finger : REMD simulation
 &cntrl
  imin = 0, ipol = 1, irest = 1, ntx = 7,
  ntb = 1,
  cut = 8, ntr=0,
  ntc = 2, ntf = 2,
  temp0 = 283.90,
  ntt = 1, tautp = 0.5,
  nstlim = 40, dt = 0.001,numexchg=10,
  ntpr = 1, ntwx = 20, ntwr = 20
 /
 &ewald
  indmeth = 3,
  diptau = 9.99,
  dipmass = 0.55
 /

My output data is as following,
 
 NSTEP = 1 TIME(PS) = 550.041 TEMP(K) = 285.47 PRESS = 0.0
 Etot = -20097.5755 EKtot = 4214.2849 EPtot = -24311.8604
 BOND = 91.2338 ANGLE = 243.6370 DIHED = 226.3469
 1-4 NB = 90.4093 1-4 EEL = 760.0888 VDWAALS = 3369.8051
 EELEC = -29093.3813 EHBOND = 0.0000 RESTRAINT = 0.0000
 Dipole convergence: rms = 0.420E+00 temperature = 0.00
 ------------------------------------------------------------------------------

 NSTEP = 2 TIME(PS) = 550.042 TEMP(K) = 286.22 PRESS = 0.0
 Etot = -20104.8694 EKtot = 4225.4455 EPtot = -24330.3149
 BOND = 88.9608 ANGLE = 245.3921 DIHED = 226.7057
 1-4 NB = 90.4598 1-4 EEL = 758.3377 VDWAALS = 3367.0800
 EELEC = -29115.9920 EHBOND = 0.0000 RESTRAINT = 0.0000
 EPOLZ = 8.7411 E3BODY = 0.0000
 Dipole convergence: rms = 0.341E+00 temperature = 83.05
 ------------------------------------------------------------------------------

 NSTEP = 3 TIME(PS) = 550.043 TEMP(K) = 287.20 PRESS = 0.0
 Etot = -20098.8096 EKtot = 4239.8778 EPtot = -24338.6874
 BOND = 85.6602 ANGLE = 245.1373 DIHED = 227.2315
 1-4 NB = 90.3593 1-4 EEL = 756.4102 VDWAALS = 3363.8731
 EELEC = -29129.7588 EHBOND = 0.0000 RESTRAINT = 0.0000
 EPOLZ = 22.3996 E3BODY = 0.0000
 Dipole convergence: rms = 0.232E+00 temperature = 137.27
 ------------------------------------------------------------------------------

...............

Thanks,

Wenfei Li

Institute of Biophysics,
Nanjing University
P. R. China

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