Problems with relaxation using NPT

Hello,

I’m focusing on a system containing 1 polymer chain and around 200 solvent molecules.

I wish to relax the system to 300 K and 1 atm using NPT with OPLS force field.

In order to accelerate the simulation I use pppm instead of ewald. The initial timestep = 1 fs

The cubic simulation box begins to shrink and its volume keeps decreasing until the “cannot compute pppm” error occured.

I then continue the relaxation with timestep = 0.1 fs, however the volume begins to increase after restarting and seems to stabilize… This is extremely strange because I only modify the timestep and the system refuses to continue the shrinkage after restarting… I tried using both the binary restart file and text lammps data file for restarting the simulation. I also tried to continue with ewald instead of pppm… But in all cases the system doesn’t retrieve the trend of volume shrinkage.

Could anyone give me some advices on the relaxation ?

Please find below my input script:

----------------------------------

include “system.in.init”

read_data “system.data”

include “system.in.settings”

variable max_it equal 100000

min_style cg # conjugate gradient (CG) algorithm
minimize 1e-25 1e-25 {max_it} {max_it} # only to first initialize simulation, turned off when continuing from restart file

dump 1 all dcd 1000 system.dcd

thermo_style custom vol density temp press # specify thermodynamic variables that you want to output
thermo 1000

timestep 1

velocity all create 300.0 4928459 rot yes dist gaussian # only to first initialize simulation, turned of when continuing from restart file

fix 1 all npt temp 300.0 300.0 100.0 iso 1.0 1.0 100.0

restart 100000 crosslinker_solvent.restart

run 500000

write_data relax_restart.data

#---------------------------------------------------

units           real
atom_style	    full
kspace_style    pppm 1.0e-4
neigh_modify    every 1 delay 0 check yes
pair_style      lj/cut/coul/long  10.0                                                                                                                     
bond_style      harmonic
angle_style     harmonic
dihedral_style  opls

There is not enough information here to make any specific statements, but it looks like 1fs is like too large a time step. Try with 0.25fs instead, and also first run for some time with nvt before switching to npt.

If you still have problems after applying those two suggestions, please provide all files and also the log file.

Thank you so much for your reply ! I’m currently testing your suggestions.

In fact as a new user I can’t upload files so I just attached the log file here for a typical cannot compute pppm error

Last night I managed to circumvent this issue with a … not so elegant method…

Each time I encountered the “cannot compute pppm error”, I relaunched the simulation with the latest restart file/data file. All the parameters are the same, including the timestep = 1 fs.

When the volume shrinked to a point where the same error occurred again, I repeated the same operations.

After 3 or 4 relaunchs, we reached the equilibrium density…

This is a pretty violent approach but it works, very efficiently. And it’s a lot faster than my previous strategy (switching to smaller timesteps or ewald). We just need to relaunch manually the same script several times…

I also saw somewhere on the forum that fix npt needs to be reinitiated if the box size changed to much, so the relaunch seems to be a way for “reinitiation”… The cause of this problem might be due to the volume shrinkage with NPT. After each relaunch the volume decreases for some time and another out of range atoms, cannot compute pppm error occurs, which means that we need to reinitiate again.

Would you mind telling me if there is another way to reinitiate regularly in case that the volume changes a lot with NPT ? So that we can avoid relaunching manually the simulation…

Thank you very much in advance.

Tianchi

#---------------------------------------------------

LAMMPS (29 Oct 2020)
using 1 OpenMP thread(s) per MPI task
Reading data file …
orthogonal box = (-69.460211 -69.460211 -69.460211) to (69.460211 69.460211 69.460211)
3 by 4 by 4 MPI processor grid
reading atoms …
6085 atoms
reading velocities …
6085 velocities
scanning bonds …
4 = max bonds/atom
scanning angles …
6 = max angles/atom
scanning dihedrals …
30 = max dihedrals/atom
reading bonds …
5868 bonds
reading angles …
11060 angles
reading dihedrals …
13548 dihedrals
Finding 1-2 1-3 1-4 neighbors …
special bond factors lj: 0.0 0.0 0.0
special bond factors coul: 0.0 0.0 0.0
4 = max # of 1-2 neighbors
10 = max # of 1-3 neighbors
24 = max # of 1-4 neighbors
26 = max # of special neighbors
special bonds CPU = 0.002 seconds
read_data CPU = 0.067 seconds
PPPM initialization …
WARNING: System is not charge neutral, net charge = -0.1123 (src/kspace.cpp:313)
using 12-bit tables for long-range coulomb (src/kspace.cpp:328)
G vector (1/distance) = 0.22178144
grid = 45 45 45
stencil order = 5
estimated absolute RMS force accuracy = 0.035979295
estimated relative force accuracy = 0.00010835058
using double precision FFTW3
3d grid and FFT values/proc = 5780 2160
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 = 24 24 24
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/newton
bin: standard
Setting up Verlet run …
Unit style : real
Current step : 0
Time step : 1.0
Per MPI rank memory allocation (min/avg/max) = 19.87 | 19.95 | 20.39 Mbytes
Volume Density Temp Press
2681009.1 0.022902972 300.52065 -3.1611842
2740055.2 0.022409431 302.65696 -4.6853735
2806570 0.021878335 297.34069 -7.1047648
3032688.9 0.020247074 298.81313 -10.044813
3176511.5 0.01933035 303.42791 -3.3296846
3055540.1 0.020095654 300.6911 -0.058853951
2997326.2 0.020485951 296.407 14.075634
3274083.4 0.01875428 298.55679 13.746628
3735347.4 0.016438385 302.81084 -3.8863322
3812437.9 0.016105987 296.46992 -1.8150647
3667081.7 0.016744399 304.74565 6.3508647
3484537.9 0.017621584 298.61493 13.795785
3346731.2 0.018347179 294.7785 20.0333
3120925.6 0.019674637 300.77613 10.471741
3102675.6 0.019790363 302.05439 -35.232131
3248502.9 0.018901962 302.30676 -23.377178
3364194.2 0.018251942 296.27398 -5.1134209
3608842.5 0.017014618 304.58194 11.977185
3744512.1 0.016398152 299.84055 13.253365
3887974.1 0.015793078 300.07742 -0.10763535
3521913.7 0.017434577 300.11982 -9.9904621
3201147.1 0.019181586 296.18828 -3.280912
2992985.1 0.020515664 294.79584 13.118351
2943881.6 0.020857862 298.57701 16.954071
3313095.1 0.018533448 301.49007 -6.5492646
3578132.6 0.017160649 298.96479 -9.2857143
3657140.5 0.016789915 303.59706 -8.8200485
3542869.6 0.017331453 302.39226 8.8521241
3396642.8 0.018077578 304.93756 17.61928
3033723.6 0.020240169 301.86696 6.4687138
3048827.6 0.020139898 298.70584 -2.9764166
3133125.5 0.019598027 296.47086 -20.864485
3199593.2 0.019190901 299.11671 -7.9055835
3337267.6 0.018399207 303.12286 11.710259
3383146.7 0.018149694 301.86764 -4.1682049
3144431.1 0.019527563 297.11776 -6.99579
3089978.9 0.019871682 302.47332 -19.00658
2975612.5 0.020635441 302.04762 6.6790717
3019720.1 0.020334029 298.08213 2.0385365
2909010.7 0.021107889 304.02623 18.942936
2729736.7 0.022494139 301.02629 10.696368
2773001.9 0.022143179 297.51814 -24.930288
3108425.9 0.019753753 294.88447 -17.076694
3171664.5 0.01935989 301.10198 7.0844147
3231364.4 0.019002214 299.8964 22.147302
3101419.3 0.01979838 295.13579 4.4185095
3096817 0.019827803 298.93985 19.410343
2983517.8 0.020580765 300.99885 8.8488713
2668762.2 0.023008073 300.70133 -2.5356531
2738303.6 0.022423765 301.43677 13.013271
2660358.2 0.023080756 297.32118 -0.90976771
2733310.3 0.02246473 299.79182 13.038183
2691781.3 0.022811317 303.19774 -19.663281
2739317 0.022415469 299.44124 1.8273545
2779143.7 0.022094243 295.56678 24.209018
2828326.3 0.02171004 299.03667 27.73251
2547024.1 0.024107772 302.89702 12.059627
2511013.1 0.024453507 300.88788 -4.5547439
2402234.7 0.025560815 299.45558 -4.9699648
2302876.7 0.02666364 295.08626 6.2095025
2031927.1 0.030219134 306.60421 13.973624
2249070.2 0.027301539 306.66525 -12.812358
2444158.2 0.025122383 299.39514 -9.8386231
2638120.1 0.023275316 300.41295 2.9094824
2732739 0.022469426 298.55295 -12.678648
2573104.7 0.02386342 295.81073 -9.4949999
2426786.2 0.025302219 299.45103 0.83562982
2361097.8 0.026006155 299.575 -10.985484
2345208.9 0.026182348 298.90689 -5.9180272
2377978.3 0.025821546 299.34903 -19.836887
2441868.2 0.025145943 298.18064 13.038139
2237944.6 0.027437265 301.4739 6.0249433
2131282.3 0.028810391 304.77243 3.2393257
2007699.8 0.030583795 301.25776 -13.261294
1873370.6 0.032776791 303.45724 9.554189
2011845 0.030520779 299.79078 20.311097
2079947.5 0.029521455 301.5644 19.093567
2309720.8 0.026584632 296.11132 10.616322
2373479.2 0.025870493 296.20942 24.703315
2366487 0.025946932 305.75052 23.044931
2160761 0.028417339 297.88198 -10.492233
1976462 0.031067169 301.62117 -9.6620792
2030444.7 0.030241197 297.17788 -14.564784
1994624.5 0.030784279 295.98805 1.0412004
1964759.7 0.031252207 300.69672 16.245966
1792036.2 0.034264418 302.98557 30.497875
1716586.3 0.035770458 296.89382 -3.3152943
1694707.3 0.03623226 296.2026 -44.017768
1682974.6 0.036484851 302.37925 -11.150764
1683420 0.036475197 300.80936 12.621726
1623328.3 0.037825421 297.17363 54.596977
1545278.6 0.039735926 300.90742 41.554463
1589673.6 0.038626216 297.3958 -50.602411
1612338.2 0.038083249 295.51277 -43.51423
1563068.5 0.039283677 300.43008 -17.103571
1516967.1 0.040477527 302.40327 33.327419
1391638 0.04412288 300.97766 -6.0814536
1287497 0.047691822 303.34015 -52.378093
1314169.3 0.046723872 299.37123 -9.3338833
1289584.3 0.047614629 297.95502 -1.367003
1173826.5 0.052310179 303.53085 51.818339
1138497.6 0.053933429 297.48112 27.86496
1137991.8 0.0539574 295.72051 -29.824657
1102813.2 0.055678585 296.90208 -38.547117
1129754.2 0.054350827 296.22987 -10.928439
1003953.2 0.061161292 304.05249 -23.786961
922179.91 0.066584704 293.52457 70.641001
841270.58 0.072988499 295.70492 3.5147258
808710.73 0.07592712 297.54007 -95.070232
758028.08 0.081003697 298.95563 -54.672823
689544.14 0.089048798 302.62801 13.143828
631819.48 0.097184527 299.2462 11.657463
610519.63 0.10057511 298.36057 -155.38154
555937.08 0.11044969 297.75804 -129.51686
526801.3 0.11655832 300.2597 73.584626
519165.89 0.11827256 295.6818 119.52518
518808.88 0.11835394 300.15692 -101.64002
494995.31 0.1240478 304.83785 -109.54322
440849.13 0.13928365 292.29278 8.8120407
416770.13 0.1473308 298.31963 64.428837
388438.91 0.15807654 302.45374 76.992332
363281.69 0.16902332 304.17474 48.135254
311745.35 0.19696549 302.33354 11.177268
293697.89 0.20906884 302.31991 -303.51287

MPI_ABORT was invoked on rank 28 in communicator MPI_COMM_WORLD
with errorcode 1.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.

ERROR on proc 28: Out of range atoms - cannot compute PPPM (src/KSPACE/pppm.cpp:1909)

You can put files into dropbox, onedrive, google drive or equivalent and provide a link to those.

Why do you think a timestep of 1fs is suitable? Your input looks like a typical molecular system and I don’t see any indication of constraining hydrogen atoms. Under those circumstances 1 fs would be too much.

Please see the manual: fix nvt command — LAMMPS documentation

As I already mentioned, you should start without the changing volume. You have a rather small system with a small number of atoms, there is is much more likely to reach unstable time integration of the nose-hoover barostat chains due to the large pressure fluctuations, especially when your initial structure is rather far away from equilibrium (and thus your pressure may not be stablilized yet).

Please also note that stress/pressure is far more sensitive to Kspace convergence than forces, so using a convergence or 1e-4 for PPPM or Ewald is not likely sufficient to get converged correct pressure. In my personal experience you may need something more like 1e-6.