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.scripps.edu
To unsubscribe, send "unsubscribe amber" to majordomo.scripps.edu
Received on Thu Sep 01 2005 - 02:53:01 PDT
