Using SPCFW, I have only been able to run a 'production run' with an
integration time step of 2fs. Lower times (1fs or .25fs) result in wild
fluctuations of temperature and energy. Here is what I have tried:
After 100ps of adjusting my density using this input file:
&cntrl
imin = 0, irest = 1, ntx = 5,
ntb = 2, pres0 = 1.0, ntp = 1,
taup = 2.0,
cut = 10, ntr = 0,
ntc = 1, ntf = 1, jfastw=4,
tempi = 300.0, temp0 = 300.0,
ntt = 3, gamma_ln = 1.0,
nstlim = 50000, dt = 0.002,
ntpr = 50, ntwx = 1000, ntwr = 1000
&end
this is the output from my last few steps:
NSTEP = 49850 TIME(PS) = 115.700 TEMP(K) = 300.96 PRESS =
83.6
Etot = -25025.9679 EKtot = 10900.7514 EPtot =
-35926.7192
BOND = 6114.6090 ANGLE = 3248.5105 DIHED =
1119.9839
1-4 NB = 373.9743 1-4 EEL = 5938.9232 VDWAALS =
7375.5304
EELEC = -60098.2505 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3183.7469 VIRIAL = 2969.3649 VOLUME =
118771.0214
Density =
1.0445
Ewald error estimate: 0.6706E-04
------------------------------------------------------------------------------
NSTEP = 49900 TIME(PS) = 115.800 TEMP(K) = 298.37 PRESS =
444.9
Etot = -25012.9862 EKtot = 10806.8796 EPtot =
-35819.8659
BOND = 5980.8253 ANGLE = 3239.7259 DIHED =
1108.0627
1-4 NB = 400.8266 1-4 EEL = 5937.7246 VDWAALS =
7399.0622
EELEC = -59886.0931 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3134.1026 VIRIAL = 1993.1803 VOLUME =
118784.5104
Density =
1.0444
Ewald error estimate: 0.1066E-03
------------------------------------------------------------------------------
NSTEP = 49950 TIME(PS) = 115.900 TEMP(K) = 298.15 PRESS =
-202.4
Etot = -25001.2443 EKtot = 10798.8065 EPtot =
-35800.0508
BOND = 6203.6212 ANGLE = 3231.7346 DIHED =
1104.1238
1-4 NB = 383.4923 1-4 EEL = 5950.3001 VDWAALS =
7354.6589
EELEC = -60027.9817 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3161.4493 VIRIAL = 3680.5986 VOLUME =
118818.7381
Density =
1.0441
Ewald error estimate: 0.1072E-03
------------------------------------------------------------------------------
NSTEP = 50000 TIME(PS) = 116.000 TEMP(K) = 299.97 PRESS =
168.5
Etot = -24972.9685 EKtot = 10864.7098 EPtot =
-35837.6783
BOND = 6084.2312 ANGLE = 3297.5199 DIHED =
1129.0296
1-4 NB = 368.0507 1-4 EEL = 5915.6890 VDWAALS =
7390.8470
EELEC = -60023.0457 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3135.4620 VIRIAL = 2703.2618 VOLUME =
118772.7244
Density =
1.0445
Ewald error estimate: 0.2171E-04
------------------------------------------------------------------------------
A V E R A G E S O V E R 50000 S T E P S
NSTEP = 50000 TIME(PS) = 116.000 TEMP(K) = 300.30 PRESS =
-147.5
Etot = -24768.9263 EKtot = 10876.6253 EPtot =
-35645.5516
BOND = 6073.7988 ANGLE = 3234.2054 DIHED =
1111.4209
1-4 NB = 382.0535 1-4 EEL = 5933.0745 VDWAALS =
7320.4364
EELEC = -59700.5411 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3158.6254 VIRIAL = 3607.4813 VOLUME =
127782.9686
Density =
0.9811
Ewald error estimate: 0.7549E-04
------------------------------------------------------------------------------
R M S F L U C T U A T I O N S
NSTEP = 50000 TIME(PS) = 116.000 TEMP(K) = 2.14 PRESS =
292.5
Etot = 438.0896 EKtot = 77.6314 EPtot =
431.1896
BOND = 95.5984 ANGLE = 53.8788 DIHED =
12.9521
1-4 NB = 8.5207 1-4 EEL = 20.6381 VDWAALS =
117.8675
EELEC = 541.8019 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 42.7825 VIRIAL = 837.2441 VOLUME =
13892.2384
Density =
0.0946
Ewald error estimate: 0.5565E-04
*****************************************
******************************************
******************************************
Then I tried running a 10ps test production run using Berendsen algorithm
and 1fs integration time but my energy and temperature fluctuated from 300K
to ~330K. I thought my time step was too large so I tried .25fs integration
time:
&cntrl
imin = 0, irest = 1, ntx = 5,
ntb = 2, pres0=1.0, ntp=1,
taup= 11.0,
cut = 10, ntr = 0,
ntc = 1, ntf = 1, jfastw=4, tol=.000001,
tempi = 300.0, temp0 = 300.0, tautp=11.0,
ntt = 1, nscm=0,
nstlim = 40000, dt = 0.00025,
ntpr = 1, ntwx = 6, ntwr = 1000,
&end
&ewald
ew_type=0, dsum_tol=0.000001
&end
This time the fluctuations went from 300K to 400K and back. Here is the
beginning of my simulation and the end of 10ps run:
NSTEP = 7 TIME(PS) = 116.002 TEMP(K) = 279.45 PRESS =
97.9
Etot = -23364.9284 EKtot = 10121.6123 EPtot =
-33486.5406
BOND = 8347.2725 ANGLE = 3507.6822 DIHED =
1129.2365
1-4 NB = 371.4838 1-4 EEL = 5926.1573 VDWAALS =
7412.6521
EELEC = -60181.0249 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3117.2778 VIRIAL = 2866.2173 VOLUME =
118773.6598
Density =
1.0445
Ewald error estimate: 0.3033E-04
------------------------------------------------------------------------------
NSTEP = 8 TIME(PS) = 116.002 TEMP(K) = 292.69 PRESS =
122.7
Etot = -23368.7147 EKtot = 10601.2506 EPtot =
-33969.9653
BOND = 7814.4207 ANGLE = 3530.2662 DIHED =
1129.9560
1-4 NB = 371.9707 1-4 EEL = 5926.9408 VDWAALS =
7415.1497
EELEC = -60158.6694 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3115.3400 VIRIAL = 2800.6133 VOLUME =
118773.6715
Density =
1.0445
Ewald error estimate: 0.3138E-04
------------------------------------------------------------------------------
NSTEP = 9 TIME(PS) = 116.002 TEMP(K) = 311.15 PRESS =
153.1
Etot = -23373.8170 EKtot = 11269.6860 EPtot =
-34643.5030
BOND = 7090.0578 ANGLE = 3547.2938 DIHED =
1130.6773
1-4 NB = 372.4193 1-4 EEL = 5927.7072 VDWAALS =
7417.4162
EELEC = -60129.0747 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3113.4464 VIRIAL = 2720.8810 VOLUME =
118773.6861
Density =
1.0445
Ewald error estimate: 0.3215E-04
------------------------------------------------------------------------------
NSTEP = 10 TIME(PS) = 116.002 TEMP(K) = 332.60 PRESS =
188.0
Etot = -23379.6520 EKtot = 12046.4421 EPtot =
-35426.0941
BOND = 6256.7042 ANGLE = 3558.4580 DIHED =
1131.3964
1-4 NB = 372.8263 1-4 EEL = 5928.4485 VDWAALS =
7419.4784
EELEC = -60093.4060 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3111.5882 VIRIAL = 2629.5405 VOLUME =
118773.7044
Density =
1.0445
Ewald error estimate: 0.3484E-04
------------------------------------------------------------------------------
NSTEP = 11 TIME(PS) = 116.003 TEMP(K) = 354.50 PRESS =
226.1
Etot = -23385.5712 EKtot = 12839.7070 EPtot =
-36225.2782
BOND = 5408.4531 ANGLE = 3563.5753 DIHED =
1132.1101
1-4 NB = 373.1899 1-4 EEL = 5929.1560 VDWAALS =
7421.3749
EELEC = -60053.1374 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3109.7569 VIRIAL = 2529.8264 VOLUME =
118773.7269
Density =
1.0445
Ewald error estimate: 0.3489E-04
------------------------------------------------------------------------------
NSTEP = 12 TIME(PS) = 116.003 TEMP(K) = 374.30 PRESS =
266.3
Etot = -23390.8793 EKtot = 13556.8356 EPtot =
-36947.7149
BOND = 4640.2589 ANGLE = 3562.5857 DIHED =
1132.8153
1-4 NB = 373.5096 1-4 EEL = 5929.8215 VDWAALS =
7423.1720
EELEC = -60009.8778 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3107.9444 VIRIAL = 2425.0798 VOLUME =
118773.7540
Density =
1.0445
Ewald error estimate: 0.3344E-04
------------------------------------------------------------------------------
NSTEP = 13 TIME(PS) = 116.003 TEMP(K) = 389.70 PRESS =
307.1
Etot = -23394.9910 EKtot = 14114.9133 EPtot =
-37509.9043
BOND = 4037.1644 ANGLE = 3555.5494 DIHED =
1133.5094
1-4 NB = 373.7864 1-4 EEL = 5930.4371 VDWAALS =
7424.8971
EELEC = -59965.2481 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3106.1441 VIRIAL = 2318.5382 VOLUME =
118773.7860
Density =
1.0445
Ewald error estimate: 0.3093E-04
------------------------------------------------------------------------------
NSTEP = 14 TIME(PS) = 116.003 TEMP(K) = 398.96 PRESS =
347.6
Etot = -23397.4239 EKtot = 14450.1729 EPtot =
-37847.5968
BOND = 3664.6946 ANGLE = 3542.6416 DIHED =
1134.1902
1-4 NB = 374.0221 1-4 EEL = 5930.9956 VDWAALS =
7426.5973
EELEC = -59920.7383 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3104.3515 VIRIAL = 2212.9321 VOLUME =
118773.8228
Density =
1.0445
Ewald error estimate: 0.3112E-04
------------------------------------------------------------------------------
NSTEP = 15 TIME(PS) = 116.004 TEMP(K) = 401.03 PRESS =
386.8
Etot = -23397.9321 EKtot = 14525.1918 EPtot =
-37923.1240
BOND = 3561.5117 ANGLE = 3524.1440 DIHED =
1134.8557
1-4 NB = 374.2189 1-4 EEL = 5931.4907 VDWAALS =
7428.2855
EELEC = -59877.6303 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3102.5650 VIRIAL = 2110.6734 VOLUME =
118773.8646
Density =
1.0445
Ewald error estimate: 0.3081E-04
------------------------------------------------------------------------------
NSTEP = 16 TIME(PS) = 116.004 TEMP(K) = 395.73 PRESS =
424.1
Etot = -23396.4697 EKtot = 14333.0527 EPtot =
-37729.5224
BOND = 3735.1627 ANGLE = 3500.4356 DIHED =
1135.5041
1-4 NB = 374.3796 1-4 EEL = 5931.9169 VDWAALS =
7429.9822
EELEC = -59836.9035 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3100.7857 VIRIAL = 2013.2774 VOLUME =
118773.9110
Density =
1.0445
Ewald error estimate: 0.2718E-04
------------------------------------------------------------------------------
NSTEP = 17 TIME(PS) = 116.004 TEMP(K) = 383.72 PRESS =
459.3
Etot = -23393.2032 EKtot = 13897.9977 EPtot =
-37291.2009
BOND = 4161.3986 ANGLE = 3471.9828 DIHED =
1136.1339
1-4 NB = 374.5070 1-4 EEL = 5932.2703 VDWAALS =
7431.6737
EELEC = -59799.1671 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3099.0182 VIRIAL = 1921.2038 VOLUME =
118773.9619
Density =
1.0445
Ewald error estimate: 0.2929E-04
------------------------------------------------------------------------------
**************************************
**************************************
************************************** and end:
NSTEP = 39999 TIME(PS) = 126.000 TEMP(K) = 333.76 PRESS =
-438.8
Etot = -24808.0715 EKtot = 12088.4960 EPtot =
-36896.5675
BOND = 5499.0118 ANGLE = 3279.2645 DIHED =
1115.4754
1-4 NB = 380.0574 1-4 EEL = 5900.1106 VDWAALS =
7426.4185
EELEC = -60496.9057 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3071.5697 VIRIAL = 4192.3129 VOLUME =
118298.8735
Density =
1.0487
Ewald error estimate: 0.3305E-04
------------------------------------------------------------------------------
NSTEP = 40000 TIME(PS) = 126.000 TEMP(K) = 334.10 PRESS =
-442.1
Etot = -24808.1869 EKtot = 12100.8138 EPtot =
-36909.0007
BOND = 5478.7897 ANGLE = 3279.1340 DIHED =
1115.5362
1-4 NB = 380.2458 1-4 EEL = 5900.4396 VDWAALS =
7425.8974
EELEC = -60489.0434 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3072.7838 VIRIAL = 4202.0443 VOLUME =
118298.8208
Density =
1.0487
Ewald error estimate: 0.3505E-04
------------------------------------------------------------------------------
A V E R A G E S O V E R 40000 S T E P S
NSTEP = 40000 TIME(PS) = 126.000 TEMP(K) = 342.99 PRESS =
-97.8
Etot = -24140.8552 EKtot = 12422.8359 EPtot =
-36563.6910
BOND = 5699.3751 ANGLE = 3192.6346 DIHED =
1113.5491
1-4 NB = 380.5053 1-4 EEL = 5946.0372 VDWAALS =
7420.4404
EELEC = -60316.2328 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3130.4151 VIRIAL = 3380.7208 VOLUME =
118484.6181
Density =
1.0470
Ewald error estimate: 0.2593E-04
------------------------------------------------------------------------------
R M S F L U C T U A T I O N S
NSTEP = 40000 TIME(PS) = 126.000 TEMP(K) = 6.30 PRESS =
285.4
Etot = 413.0053 EKtot = 228.3143 EPtot =
285.6388
BOND = 190.4687 ANGLE = 59.5774 DIHED =
12.1612
1-4 NB = 8.0443 1-4 EEL = 20.2781 VDWAALS =
119.5423
EELEC = 240.5751 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 43.3306 VIRIAL = 734.2630 VOLUME =
143.7660
Density =
0.0013
Ewald error estimate: 0.1137E-04
------------------------------------------------------------------------------
****************************************
*****************************************
***************************************
I tried each of the following separately: I got ride of tol and dsum_tol
values, changed nscm to the default, changed taup and tautp to default, I
tried Langevin algorithm with ln-gamma =1, and used NVE ensemble. NVE
reduced the oscillations to 330K and the entire system seems to oscillate
like a dampened spring. By the end of the simulation the oscillation is from
305K to 310K.
The only thing that worked is changing the integration time to 2fs. This is
the input file I used:
&cntrl
imin = 0, irest = 1, ntx = 5,
ntb = 2, pres0=1.0, ntp=1,
taup= 11.0,
cut = 10, ntr = 0,
ntc = 1, ntf = 1, jfastw=4, tol=.000001,
tempi = 300.0, temp0 = 300.0, tautp=11.0,
ntt = 1, nscm=0,
nstlim = 5000, dt = 0.002,
ntpr = 1, ntwx = 1, ntwr = 1000,
&end
&ewald
ew_type=0, dsum_tol=0.000001
&end
The output still oscillates but its between 300K and ~303K. For example:
NSTEP = 44 TIME(PS) = 116.088 TEMP(K) = 298.42 PRESS =
-166.5
Etot = -24967.9730 EKtot = 10808.6641 EPtot =
-35776.6371
BOND = 6213.4088 ANGLE = 3293.1923 DIHED =
1103.4504
1-4 NB = 362.7627 1-4 EEL = 5927.8977 VDWAALS =
7306.0860
EELEC = -59983.4351 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3100.1324 VIRIAL = 3527.1050 VOLUME =
118765.2791
Density =
1.0446
Ewald error estimate: 0.3055E-04
------------------------------------------------------------------------------
NSTEP = 45 TIME(PS) = 116.090 TEMP(K) = 299.19 PRESS =
-275.5
Etot = -24996.9683 EKtot = 10836.6822 EPtot =
-35833.6505
BOND = 6201.4941 ANGLE = 3294.0531 DIHED =
1104.3435
1-4 NB = 363.7658 1-4 EEL = 5932.0194 VDWAALS =
7301.6649
EELEC = -60030.9913 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3104.6453 VIRIAL = 3811.1847 VOLUME =
118765.1178
Density =
1.0446
Ewald error estimate: 0.2801E-04
------------------------------------------------------------------------------
NSTEP = 46 TIME(PS) = 116.092 TEMP(K) = 299.74 PRESS =
-210.9
Etot = -24964.1320 EKtot = 10856.3657 EPtot =
-35820.4977
BOND = 6119.1690 ANGLE = 3253.6060 DIHED =
1105.4888
1-4 NB = 366.7742 1-4 EEL = 5936.0300 VDWAALS =
7291.0675
EELEC = -59892.6332 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3110.7166 VIRIAL = 3651.4493 VOLUME =
118764.8515
Density =
1.0446
Ewald error estimate: 0.2738E-04
------------------------------------------------------------------------------
NSTEP = 47 TIME(PS) = 116.094 TEMP(K) = 303.03 PRESS =
-138.6
Etot = -25011.1851 EKtot = 10975.4574 EPtot =
-35986.6425
BOND = 5923.2753 ANGLE = 3223.4300 DIHED =
1106.5644
1-4 NB = 369.9693 1-4 EEL = 5938.1940 VDWAALS =
7276.9319
EELEC = -59825.0076 EHBOND = 0.0000 RESTRAINT =
0.0000
EKCMT = 3117.2980 VIRIAL = 3472.7293 VOLUME =
118764.6474
Density =
1.0446
Ewald error estimate: 0.3496E-04
------------------------------------------------------------------------------
******************************************
*****************************************
******************************************
The original paper by Wu, Tepper and Voth (J. Chem. Phys. 124, 024503
(2006)) mentions that they used NPT at 298.16K, 1 atm, Nose-Hoover
thermostat and barostat, PME of 10^-6, cut of 9 angstroms and integration
step of 1fs. They performed their simulation in DL_POLY 2.13.
I am not sure what I am doing wrong. Any help would be appreciated.
Thanks
Naser
-----------------------------------------------------------------------
The AMBER Mail Reflector
To post, send mail to amber.scripps.edu
To unsubscribe, send "unsubscribe amber" (in the *body* of the email)
to majordomo.scripps.edu
Received on Fri Dec 05 2008 - 15:49:05 PST
