Hello,
I'm encountering severe problems with the tip4p potential for water molecules. Even a simple relaxation does not seem to work, producing some distorted water molecules in which one H apparently switched with O. Moreover, total energy and temperature grow during nve dynamics. Input files are attached. Tested with lammps-12Dec18.
I will file a bug report in case no one sees problems here...
Thank you,
Roberto
water.dump (419 KB)
relaxwater.in (1.43 KB)
Hello,
I'm encountering severe problems with the tip4p potential for water molecules. Even a simple relaxation does not seem to work, producing some distorted water molecules in which one H apparently switched with O. Moreover, total energy and temperature grow during nve dynamics. Input files are attached. Tested with lammps-12Dec18.
I will file a bug report in case no one sees problems here...
i can run those inputs just fine and there is nothing that seems off.
visualization looks ok in VMD, after i added a suitable dump command
and wrote out a trajectory with image flags and converted the data
file into a servicable .psf file.
don't think there is a bug in any of the functionality you are using.
we would have detected it in our regression testing.
axel.
Dear Axel,
could you send me the output of your test? I've spent two days trying to figure out why is not working in my case.
I attach the xyz file containing the two frames saved by the script. You can see from the attached VMD screenshot that the second frame (relaxation with tip4p) is messed up.
Thank you,
Roberto
01+02.xyz.tar.bz2 (220 KB)
Here is the output of lammps with further 10000 nve steps showing temperature increase
lammps.out (138 KB)
this is a visualization problem and due to the fact, that you are
visualizing the system without using the image flag information.
below is the first chunk of my log file. i am also attaching a
snapshot from a VMD visualization from loading the topology first and
then a dump with image flags enabled.
in VMD then i ran the commands:
pbc box
pbc wrap -compound fragment
and created the image.
axel.
LAMMPS (4 Jan 2019)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread.
(src/comm.cpp:87)
using 1 OpenMP thread(s) per MPI task
units metal # eV, Ang, ps
variable kcal_mol equal 43.36e-3 # eV
atom_style full
newton on
boundary p p p
bond_style harmonic
angle_style harmonic
dihedral_style none
improper_style none
read_data water.dump
orthogonal box = (0.02645 0.02645 0.02641) to (35.5328 35.5328 35.4736)
1 by 1 by 1 MPI processor grid
reading atoms ...
4500 atoms
scanning bonds ...
2 = max bonds/atom
scanning angles ...
1 = max angles/atom
reading bonds ...
3000 bonds
reading angles ...
1500 angles
2 = max # of 1-2 neighbors
1 = max # of 1-3 neighbors
1 = max # of 1-4 neighbors
2 = max # of special neighbors
velocity all set 0.0 0.0 0.0
group water type 1 2
4500 atoms in group water
# TIP4P/Ice
bond_coeff 1 44.0 0.9572 # 1000 kcal/mol ~= 44 eV
angle_coeff 1 44.0 104.52
set type 1 charge -1.1794
1500 settings made for charge
set type 2 charge 0.5897
3000 settings made for charge
thermo_style custom step temp epair etotal ndanger
thermo 10
neighbor 4.0 bin
neigh_modify every 1 delay 0 check yes
timestep 1e-3 # ps
##### MINIMIZING BOND LENGTH AND ANGLE TO EQUILIBRIUM ONES ######
pair_style lj/cut 10.0
pair_coeff * * 0.00000 0.0000
min_style cg
minimize 1e-15 1e-15 100000 10000000
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 14
ghost atom cutoff = 14
binsize = 7, bins = 6 6 6
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d/newton
bin: standard
Per MPI rank memory allocation (min/avg/max) = 25.68 | 25.68 | 25.68 Mbytes
Step Temp E_pair TotEng Ndanger
0 0 0 734.12446 0
10 0 0 0.11575358 0
20 0 0 0.11574452 0
30 0 0 0.00032672595 0
40 0 0 9.9404306e-15 0
45 0 0 1.6733136e-21 0
Loop time of 12.5648 on 1 procs for 45 steps with 4500 atoms
98.5% CPU use with 1 MPI tasks x 1 OpenMP threads
Minimization stats:
Stopping criterion = linesearch alpha is zero
Energy initial, next-to-last, final =
734.124460728 1.67331364543e-21 1.67331364543e-21
Force two-norm initial, final = 559.975 9.48162e-10
Force max component initial, final = 8.47879 3.21812e-11
Final line search alpha, max atom move = 0.5 1.60906e-11
Iterations, force evaluations = 45 186
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
this is a visualization problem and due to the fact, that you are
visualizing the system without using the image flag information.
you are right, that's a visualization problem
below is the first chunk of my log file
I can't see much from it: is temperature going below 0.1 in the next?
Below is the output of nve dynamics *after relaxation*. Why T grows?
Step Temp E_pair TotEng Ndanger
10200 1.7235414 -1239.1801 -1238.5119 0
10300 3.4059539 -1345.3689 -1344.0486 0
10400 6.0789104 -1339.4616 -1337.1051 0
10500 10.140101 -1236.8799 -1232.9491 0
10600 15.000399 -1260.1217 -1254.3068 0
10700 19.810891 -1250.9076 -1243.2279 0
10800 26.390378 -1225.8311 -1215.6008 0
10900 31.390059 -1086.3953 -1074.2269 0
11000 37.377696 -1129.389 -1114.8995 0
11100 41.514043 -940.12626 -924.03329 0
11200 46.465468 -1230.6618 -1212.6494 0
11300 50.656308 -1224.0656 -1204.4286 0
11400 56.546815 -1336.3438 -1314.4233 0
11500 62.87552 -1167.2232 -1142.8495 0
11600 64.351034 -1157.1323 -1132.1866 0
11700 72.649199 -1235.0322 -1206.8696 0
11800 74.718706 -992.54822 -963.58343 0
11900 78.891951 -1218.3901 -1187.8076 0
12000 83.394732 -1437.9534 -1405.6253 0
12100 86.770017 -1308.4685 -1274.832 0
12200 93.55789 -1253.7747 -1217.5068 0
12300 95.921963 -1084.0412 -1046.8569 0
12400 99.106223 -1176.9367 -1138.518 0
12500 104.69871 -1118.7964 -1078.2098 0
12600 105.83249 -893.52033 -852.49424 0
12700 110.97307 -872.72405 -829.70521 0
12800 112.1958 -956.86308 -913.37025 0
12900 113.10854 -746.08153 -702.23488 0
13000 117.67491 -1015.6001 -969.98324 0
13100 122.38417 -848.37582 -800.93345 0
13200 121.74574 -927.4422 -880.24732 0
13300 125.49768 -898.92126 -850.27194 0
13400 127.58029 -828.14075 -778.6841 0
13500 130.99805 -767.11478 -716.33323 0
13600 132.62781 -889.44345 -838.03012 0
13700 138.58024 -795.27914 -741.55834 0
13800 140.37074 -897.16164 -842.74676 0
13900 144.03276 -560.33487 -504.5004 0
14000 144.47243 -844.46617 -788.46126 0
14100 145.24889 -920.11659 -863.81069 0
14200 151.99748 -876.82265 -817.90065 0
14300 154.85262 -761.88135 -701.85255 0
14400 154.1265 -753.77588 -694.02856 0
14500 158.13251 -840.52308 -779.22283 0
14600 162.82706 -728.22929 -665.10919 0
14700 160.25792 -906.99251 -844.86834 0
14800 164.9289 -800.10868 -736.1738 0
14900 167.10675 -859.16171 -794.38258 0
15000 167.2073 -810.98546 -746.16735 0
15100 167.59522 -510.73948 -445.77099 0
15200 171.31108 -602.92327 -536.51433 0
15300 173.86829 -738.77332 -671.37308 0
15400 178.89237 -556.6632 -487.31536 0
15500 179.86531 -791.89961 -722.17461 0
15600 183.60266 -871.41951 -800.24572 0
15700 182.60102 -769.80071 -699.01521 0
15800 186.77408 -867.3508 -794.94761 0
15900 185.46399 -935.35008 -863.45475 0
16000 192.23635 -890.95312 -816.43247 0
16100 189.18661 -683.85016 -610.51175 0
16200 191.70509 -672.27607 -597.96137 0
16300 194.05187 -916.17798 -840.95355 0
16400 195.46799 -1087.347 -1011.5736 0
16500 198.19241 -875.98228 -799.15276 0
16600 203.37121 -894.36118 -815.5241 0
16700 203.64448 -742.04846 -663.10544 0
16800 201.7826 -856.83645 -778.61519 0
16900 206.64605 -978.32264 -898.21606 0
17000 209.79908 -1113.586 -1032.2571 0
17100 212.40683 -862.91374 -780.57399 0
17200 212.90919 -805.61419 -723.0797 0
17300 213.33858 -682.53206 -599.8311 0
17400 213.5582 -834.38365 -751.59756 0
17500 221.4332 -882.84133 -797.00249 0
17600 222.64367 -964.52552 -878.21744 0
17700 223.47474 -972.15028 -885.52004 0
17800 220.27395 -884.84325 -799.4538 0
17900 226.84502 -943.20454 -855.2678 0
18000 228.09172 -742.11101 -653.69099 0
18100 230.44844 -578.49201 -489.15841 0
18200 228.34001 -780.08807 -691.5718 0
18300 234.12934 -980.92639 -890.16588 0
18400 234.47744 -1015.8595 -924.96407 0
18500 235.25719 -944.39151 -853.19378 0
18600 240.83149 -885.72072 -792.36211 0
18700 241.42545 -637.31411 -543.72525 0
18800 239.17654 -768.96606 -676.249 0
18900 243.96214 -778.04473 -683.47252 0
19000 251.94855 -575.42154 -477.75339 0
19100 251.8109 -819.04944 -721.43464 0
19200 247.12142 -743.25662 -647.45971 0
19300 254.5733 -841.19349 -742.50785 0
19400 255.56345 -852.89153 -753.82206 0
19500 256.17926 -820.49287 -721.18468 0
19600 263.95489 -842.05036 -739.72793 0
19700 259.6449 -818.40644 -717.75479 0
19800 258.06728 -743.89087 -643.85078 0
19900 259.24343 -964.52923 -864.03321 0
20000 262.83709 -1060.9175 -959.02843 0
> this is a visualization problem and due to the fact, that you are
> visualizing the system without using the image flag information.you are right, that's a visualization problem
> below is the first chunk of my log file
I can't see much from it: is temperature going below 0.1 in the next?
Below is the output of nve dynamics *after relaxation*. Why T grows?
equilibration. note that your potential energy (E_pair) goes down at
the same time.
please also note, that you are using an unshifted and unsmoothed
cutoff coulomb interaction.
that cannot result in good energy conservation. when the system equilibrates.
total energy should be (mostly) conserved. it may be less drastic after adding
pair_modify shift yes
but ultimately, you need to use long-range coulomb to have good
(total) energy conservation after the system has equilibrated,
axel.
Dear Axel,
thanks a lot for your suggestions. I’ve tried using smoother potentials and it indeed helps.
However, with the settings below I still obtain a T increase of ~2K/ps, which seems incompatible with microcanonic simulations in the ns timescale.
Could this residual grow of T be due to the shake constraint?
Thank you
Roberto
pair_style hybrid/overlay tip4p/long 1 2 1 1 0.1577 8.5 lj/smooth/linear 10.0
pair_coeff * * lj/smooth/linear 0.00000 0.0000
pair_coeff 1 1 lj/smooth/linear $(0.21084*v_kcal_mol) 3.1668
pair_coeff * * tip4p/long
kspace_style pppm/tip4p 1e-6
fix 1 all shake 1e-4 20 500 b 1 a 1 t 1 2
fix 2 all nve
After more tests I withdraw my previous statement: I see good energy conservation and T stabilizes at some level depending on annealing and Coulomb accuracy.
Solved, thank you again.
Roberto