At the late stage of NVT calculation, the atoms no longer move

LAMMPS LAMMPS General Discussion
1

I used the system after NPT and wanted to calculate the diffusion coefficient of the atom under NVT. But when I get the results, it is found that at the beginning of the reaction, the atoms are still moving, but at the end of the process, the atoms are not moving. So I would like to ask what the reason is and how to modify it.

Here is my code:
units real
atom_style full
pair_style lj/charmm/coul/long 9.0 10.0
bond_style harmonic
angle_style harmonic
kspace_style pppm 0.0001

read_data “system_after_npt.data”

bond_coeff 1 600.0 1.0
angle_coeff 1 75.0 109.47
pair_coeff 1 1 0.1553 3.166
pair_coeff 2 2 0.0 0.0
group spce type 1 2
group na type 3
group cl type 4
pair_coeff 3 3 0.3526418 2.1595384928
pair_coeff 4 4 0.0127850 4.8304528498

fix fRattleSPCE spce rattle 0.0001 10 100 b 1 a 1

timestep 1.0
dump 1 all custom 500 traj_nvt.lammpstrj id mol type x y z ix iy iz
fix fxnvt all nvt temp 300.0 300.0 0.01 tchain 1
thermo 500

compute msdh2o spce msd
fix fxmsdh2o spce ave/time 100 5 1000 c_msdh2o[*] file h2omsd.dat

compute msdna na msd
fix fxmsdna na ave/time 100 5 1000 c_msdna[*] file namsd.dat

compute msdcl cl msd
fix fxmsdcl cl ave/time 100 5 1000 c_msdcl[*] file clmsd.dat

run 60000

write_data system_after_nvt.data

2

Hello,
I dont know anything about you system, but I am gonna assume that its salt in water at ambient temperature. If so, the damping parameter is ridiculously small, why did you choose 0.01 fs?
Best,
Simon

3

Thans for your reply, @simongravelle. And your recommendation is vary helpful. According to your suggestion, I have modified the timestep to 0.01 and get an ideal result.

But I am still confused about the result.

  1. Theoretically, the temperature of my system (Na ions, Cl ions, and water) has been maintained at 300 K, which means that the molecule should have a speed at any time.

  2. But when we extend the simulation time, neither the molecules nor the ions are moving anymore, which means that the molecules are losing kinetic energy.

Whether these two points are in conflict with each other?

In addition, I have also checked the list of temperature during the process and find that temperature fluctuated around 18 K and did not reach 300 K we set.

Below is part of the log file:

Step Temp E_pair E_mol TotEng Press
0 18.92705 -12755.387 0 -12641.705 56926475
4000 19.044732 -13213.884 0 -13099.495 41699.529
8000 17.457558 -13177.384 0 -13072.528 41545.757
10000 16.614273 -13183.145 0 -13083.354 42112.301
20000 18.636657 -13219.803 0 -13107.865 42545.346
30000 18.018571 -13262.58 0 -13154.354 41655.7
60000 18.923763 -13323.119 0 -13209.456 43366.757

4

A 1.0 fs timestep should be suitable for the kind of system you have.

The 18K is only one concern, but rather the gigantic pressure is very worrisome too.
Something is gravely wrong with your system from the very beginning.
You said that you have started from a simulation that was equilibrated with fix npt.
It looks like something went wrong there already or when creating the data file.

This is not something that can be easily fixed from remote.

As a point of reference, here is a simple input deck for a tiny system demonstrating that this kind of simulation should be working fine with a 1 fs timestep and will produce proper temperature and pressure.

The molecule file, H2O.txt, is taken from the examples/mc/ folder of the LAMMPS source distribution.

units real
atom_style full
atom_modify map array
region box block -5 5 -5 5 -5 5
boundary p p p
create_box 4 box bond/types 1 angle/types 1 &
    extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2

mass 1 15.9994
mass 2 1.008
mass 3 22.989769
mass 4 35.453

pair_style lj/cut/coul/long 10.0
pair_coeff 1 1 0.1553 3.166
pair_coeff 2 2 0.0    1.0
pair_coeff 3 3 0.3526418 2.1595384928
pair_coeff 4 4 0.0127850 4.8304528498

kspace_style pppm 0.00001

bond_style harmonic
bond_coeff * 1000.0 1.0

angle_style harmonic
angle_coeff * 100.0 109.47

molecule water H2O.txt

# fill the box with randomly placed water molecules w/o overlap
create_atoms 0 random 32 34564 NULL mol water 25367 overlap 1.33
create_atoms 3 random 2  42343 NULL overlap 1.33
create_atoms 4 random 2  42343 NULL overlap 1.33

set type 1 charge -0.8476
set type 2 charge  0.4238
set type 3 charge  1.0
set type 4 charge -1.0

velocity all create 300.0 5463576
timestep 1.0

minimize 0.0 0.0 1000 10000

reset_timestep 0

# since the trajectory and forces are kept identical through fix alchemy,
# we can do fix npt simulations, but we must use the "mixed" pressure

fix integrate all nvt temp 300 300 1.0
fix rigid all rattle 0.001 10 10000 b 1 a 1

thermo_style custom step temp press etotal density pe ke

thermo 1000
run 20000

And the corresponding output:

  using 1 OpenMP thread(s) per MPI task
Created orthogonal box = (-5 -5 -5) to (5 5 5)
  1 by 2 by 2 MPI processor grid
Read molecule template water:
  1 molecules
  0 fragments
  3 atoms with max type 2
  2 bonds with max type 1
  1 angles with max type 1
  0 dihedrals with max type 0
  0 impropers with max type 0
Created 96 atoms
  using lattice units in orthogonal box = (-5 -5 -5) to (5 5 5)
  create_atoms CPU = 0.006 seconds
Created 2 atoms
  using lattice units in orthogonal box = (-5 -5 -5) to (5 5 5)
  create_atoms CPU = 0.000 seconds
Created 2 atoms
  using lattice units in orthogonal box = (-5 -5 -5) to (5 5 5)
  create_atoms CPU = 0.000 seconds
Setting atom values ...
  32 settings made for charge
Setting atom values ...
  64 settings made for charge
Setting atom values ...
  2 settings made for charge
Setting atom values ...
  2 settings made for charge
PPPM initialization ...
  using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
  G vector (1/distance) = 0.31510913
  grid = 10 10 10
  stencil order = 5
  estimated absolute RMS force accuracy = 0.0017608851
  estimated relative force accuracy = 5.3028531e-06
  using double precision KISS FFT
  3d grid and FFT values/proc = 2448 300
Generated 6 of 6 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
  update: every = 1 steps, delay = 0 steps, check = yes
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 12
  ghost atom cutoff = 12
  binsize = 6, bins = 2 2 2
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair lj/cut/coul/long, perpetual
      attributes: half, newton on
      pair build: half/bin/newton
      stencil: half/bin/3d
      bin: standard
Setting up cg style minimization ...
  Unit style    : real
  Current step  : 0
Per MPI rank memory allocation (min/avg/max) = 8.282 | 8.285 | 8.287 Mbytes
   Step          Temp          E_pair         E_mol          TotEng         Press     
         0   300            16090.637      3.7051245e-05  16179.167      4512151.3    
       422   300           -840.19853      13.578801     -738.08967      297.17441    
Loop time of 0.383567 on 4 procs for 422 steps with 100 atoms

98.6% CPU use with 4 MPI tasks x 1 OpenMP threads

Minimization stats:
  Stopping criterion = linesearch alpha is zero
  Energy initial, next-to-last, final = 
      16090.6370375965  -826.619728010597  -826.619728010597
  Force two-norm initial, final = 157158.94 0.90721236
  Force max component initial, final = 102725.86 0.42986695
  Final line search alpha, max atom move = 1.2719297e-08 5.4676055e-09
  Iterations, force evaluations = 422 871

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 0.11493    | 0.19045    | 0.2644     |  12.2 | 49.65
Bond    | 0.00086148 | 0.00091512 | 0.0010015  |   0.0 |  0.24
Kspace  | 0.075517   | 0.14949    | 0.22512    |  13.8 | 38.97
Neigh   | 0.0019501  | 0.0019727  | 0.0019948  |   0.0 |  0.51
Comm    | 0.034833   | 0.035112   | 0.035401   |   0.1 |  9.15
Output  | 0          | 0          | 0          |   0.0 |  0.00
Modify  | 0          | 0          | 0          |   0.0 |  0.00
Other   |            | 0.005624   |            |       |  1.47

Nlocal:             25 ave          30 max          23 min
Histogram: 2 1 0 0 0 0 0 0 0 1
Nghost:         2825.5 ave        2872 max        2786 min
Histogram: 1 1 0 0 0 0 1 0 0 1
Neighs:         9053.5 ave       11609 max        5796 min
Histogram: 1 0 0 1 0 0 0 0 0 2

Total # of neighbors = 36214
Ave neighs/atom = 362.14
Ave special neighs/atom = 1.92
Neighbor list builds = 6
Dangerous builds = 0
Finding RATTLE clusters ...
       0 = # of size 2 clusters
       0 = # of size 3 clusters
       0 = # of size 4 clusters
      32 = # of frozen angles
  find clusters CPU = 0.000 seconds
PPPM initialization ...
  using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
  G vector (1/distance) = 0.31510913
  grid = 10 10 10
  stencil order = 5
  estimated absolute RMS force accuracy = 0.0017608851
  estimated relative force accuracy = 5.3028531e-06
  using double precision KISS FFT
  3d grid and FFT values/proc = 2448 300
Generated 6 of 6 mixed pair_coeff terms from geometric mixing rule
Setting up Verlet run ...
  Unit style    : real
  Current step  : 0
  Time step     : 1
RATTLE stats (type/ave/delta/count) on step 0
Bond:   1   1.01151   0.013248          64
Angle:  1   107.417   3.36471           32
Per MPI rank memory allocation (min/avg/max) = 8.407 | 8.41 | 8.412 Mbytes
   Step          Temp          Press          TotEng        Density         PotEng         KinEng    
         0   314.31925     -1639.9431     -777.42449      1.1513817     -840.19853      62.774041    
      1000   301.06544      5267.2754     -677.99582      1.1513817     -738.12289      60.127065    
      2000   336.37293     -1012.8962     -700.49593      1.1513817     -767.6744       67.178475    
      3000   286.79017     -1409.7737     -692.29576      1.1513817     -749.57185      57.276091    
      4000   315.46207      3156.1877     -698.07058      1.1513817     -761.07286      63.002278    
      5000   292.66803     -614.08452     -697.48419      1.1513817     -755.93417      58.449982    
      6000   309.44319     -1149.598      -686.7636       1.1513817     -748.56383      61.800222    
      7000   254.05414      1614.5969     -689.78701      1.1513817     -740.52525      50.738237    
      8000   310.27637     -4042.8915     -680.18928      1.1513817     -742.15589      61.966619    
      9000   286.39329     -5025.28       -704.61555      1.1513817     -761.81238      57.196827    
RATTLE stats (type/ave/delta/count) on step 10000
Bond:   1   0.999997  2.24538e-05       64
Angle:  1   109.47    0.00193684        32
     10000   296.24172     -3953.4996     -697.67943      1.1513817     -756.84313      59.1637      
     11000   323.56059     -1115.9031     -690.46975      1.1513817     -755.08941      64.619668    
     12000   284.75731      3496.3049     -694.99612      1.1513817     -751.86622      56.870098    
     13000   262.9934       1931.4308     -699.40834      1.1513817     -751.93187      52.523535    
     14000   274.22344     -1159.0547     -712.0239       1.1513817     -766.79023      54.766335    
     15000   342.64355      111.45531     -701.31229      1.1513817     -769.7431       68.430807    
     16000   334.10769     -2714.1437     -690.65991      1.1513817     -757.38599      66.726074    
     17000   321.14279      1544.3886     -692.46087      1.1513817     -756.59766      64.136799    
     18000   275.20771     -3544.9569     -708.25972      1.1513817     -763.22263      54.962907    
     19000   347.23547     -1441.2803     -677.07697      1.1513817     -746.42485      69.34788     
RATTLE stats (type/ave/delta/count) on step 20000
Bond:   1   0.999989  5.30216e-05       64
Angle:  1   109.47    0.00690826        32
     20000   261.91935      800.4894      -686.0143       1.1513817     -738.32333      52.309032    
Loop time of 7.45233 on 4 procs for 20000 steps with 100 atoms

Performance: 231.874 ns/day, 0.104 hours/ns, 2683.725 timesteps/s, 268.373 katom-step/s
98.5% CPU use with 4 MPI tasks x 1 OpenMP threads

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 2.0428     | 2.8866     | 3.5676     |  32.4 | 38.73
Bond    | 0.0028994  | 0.0031749  | 0.0038875  |   0.7 |  0.04
Kspace  | 1.7449     | 2.4415     | 3.2911     |  35.7 | 32.76
Neigh   | 0.1326     | 0.13281    | 0.13304    |   0.1 |  1.78
Comm    | 0.83983    | 0.84494    | 0.85833    |   0.8 | 11.34
Output  | 0.00046679 | 0.00049966 | 0.00057587 |   0.0 |  0.01
Modify  | 1.0706     | 1.0743     | 1.0795     |   0.3 | 14.42
Other   |            | 0.06855    |            |       |  0.92

Nlocal:             25 ave          28 max          22 min
Histogram: 1 0 0 1 0 0 1 0 0 1
Nghost:        2898.25 ave        2944 max        2847 min
Histogram: 1 0 0 1 0 0 0 0 1 1
Neighs:         9061.5 ave       11597 max        5198 min
Histogram: 1 0 0 0 0 0 1 1 0 1

Total # of neighbors = 36246
Ave neighs/atom = 362.46
Ave special neighs/atom = 1.92
Neighbor list builds = 469
Dangerous builds = 0
Total wall time: 0:00:07
5

Hi, I was not suggesting you to change the timestep, which was fine, but rather the damping time of the thermostat.
Simon

6

Thanks for you suggestion. I have tried to change the parameter of Tdamp to 0.01. But it cann’t run normally, and there is an error.

ERROR: Non-numeric pressure - simulation unstable (…/fix_nh.cpp:1059)
Last command: run 10000

And then I change it to 100. In addition, I found I didn’t the parameter of velocity. So I add velocity to the run file and then get a progressive result. The result of this simulation shows that the temperature can increase but it still cann’t maintain at 300 K.

The log file was shown as following:

   Step          Temp          Press          TotEng        Density         PotEng         KinEng    
    
    90   199.10834     -7487.5881     -13259.012      1.0268299     -14454.923      1195.9109    
   100   138.83309     -2624.5806     -12103.644      1.0264893     -12937.522      833.87768    
   200   193.20207      1900.4363     -11853.246      0.99433258    -13013.681      1160.4358    
   300   227.44918      5703.7642     -11595.812      0.95356279    -12961.947      1366.1354    
   400   242.25818      6079.9469     -11358.675      0.91377723    -12813.758      1455.0831    
   500   247.79261      6438.3445     -11140.634      0.87721254    -12628.958      1488.3247    
   600   256.45293      5676.4382     -10924.448      0.84462085    -12464.79       1540.3415    
   700   269.30281      4647.9712     -10701.655      0.81477263    -12319.177      1617.5222    
   800   286.15101      4437.9799     -10462.641      0.78709833    -12181.359      1718.718            
   900   320.19948      2889.0216     -10186.82       0.76075861    -12110.044      1923.2245         
  1000   365.38588      1575.1399     -9888.3718      0.73671841    -12083.001      2194.629        
  2000   3063.7416     -17678.554      4190.4582      0.80652513    -14211.392      18401.85            
  3000   4221.0324     -19297.891      11221.292      1.2007931     -14131.631      25352.923           
  4000   4935.3441     -20035.604      16114.617      1.5078896     -13528.699      29643.316       
  5000   5610.1425     -25995.835      20521.777      1.5984796     -13174.603      33696.38         
  6000   5931.5739     -14603.189      22627.433      1.6306925     -12999.57       35627.004       
  7000   6189.3558     -20294.021      24379.513      1.6503844     -12795.814      37175.327         
  8000   6373.3455     -13303.688      25569.651      1.6610185     -12710.78       38280.431        
  9000   6330.9712     -9464.7777      25338.888      1.6644599     -12687.028      38025.917        
 10000   6324.4942     -9784.8187      25308.584      1.6677476     -12678.429      37987.014

There are two information we can get.

  1. The temperature cann’t be maintained at 300K
  2. The density is abnormal.

I want to know what causes this phenomenon and how to solve it.

Later, I simulated NVT with the results obtained from NPT, and found that after adding Velocity parameter, the temperature could be maintained at 300K, but the atoms still did not move in the later period of simulation. So I wonder if it’s because the system is so dense. Because this sodium chloride solution has a density of 1.66.

7

Thanks for your reply and the sample you provided.

The problem of temperature has been solved because I didn’t set the initial velocity of the system.

But there is another problem for NPT process.

The density of the model is about 1.66 g/cm3 for NaCl solution, and the atoms didn’t move in the NPT process. Whether it’s so dense that the atoms can’t move?

8

The simulation will do what you tell it to do.

You are asking a lot of questions of the type “my simulation does not do what I want it to do”.
That, however, is not a LAMMPS problem, but a user problem. You will have to teach yourself how to properly debug and narrow down whether everything you do in your input is suitable for what you want to achieve. We might give some advice, if you provide a suitable, simple, and small demo input, but don’t expect too much. If this turns into a “do my work and my thinking for me” scam (and it is already on the verge of becoming that), you will see interest in helping you quickly drop to zero, even for legitimate questions.

There is no point in speculating on such issues. It is trivial to answer this for yourself by setting up multiple simulations at different densities.

Please also recall that your assessment of the situation at the beginning was incorrect. Thus chances are high that it is again.