[lammps-users] Initializing velocities

Dear LAMMPS list readers,

I am new to LAMMPS and need some guidance in setting up my atomistic simulations. I’ve searched the mail archive and tried to find the solution, but to no avail. I have previously completed test runs of my simulated water-mineral system with Materials Studio, but want to transfer it to LAMMPS in order to run longer MD calculations. I have used the msi2lmp tool to create the data file from the last frame of MatStudio trajectory, but need to establish some initial velocities before running the MD that is going to generate useful data.

I’ve tried a couple of iterations and approaches. To limit the number of inscripts to reproduce, I have re-created the inscript with #OPT N# blocks to indicate the changes from one iteration to the next.

I used the velocity command to generate initial velocities followed by a fix water shake and fix all temp/rescale for a short run (OPT 1). This executes and closes normally writing all of the data. However, as stated in the LAMMPS manual, the fix temp/rescale… doesn’t update the atom positions- the atom positions are unchanging throughout the simulation. Not surprisingly, the SHAKE stats (type/ave/delta) are constant.

If I want to update the atom positions, it is advised to use a separate time integration fix like fix nve (OPT 2). If I do include the time integration fix nve, the atom positions do change and the temperature remains at the desired level (300K). The SHAKE stat ave is constant and the delta changes a bit from step to step. But the TotEng decreases (becomes more negative) while the pressure increases (becomes more positive) and visually (via VMD) the mineral collapses and LAMMPS reports and ERROR: Dump dcd of non-matching # of atoms.

If I try to bypass the velocity-scaling and simply use fix nvt (OPT 3) to try to set the temp*,* then the temperature, TotEng, and Pressure all increase until the standard result is “nan” and eventually the run crashes. If I comment out the velocity assignment while using OPT 3 and read in the last restart file from OPT 1 (read_restart…), the result at step 1 is “nan” for the SHAKE stats and for the themo output.

I would appreciate any advice, tips, or suggestions that could be provided on how to assign the initial velocities for an MD run.

Thank you,

Aric

Inscript
units real
atom_style full

pair_style lj/cut/coul/long 13.0
bond_style harmonic
angle_style harmonic
kspace_style ewald 0.0001

# Atom Definition

read_data NaMMT_vsMD_nB.lammps05

# Settings

# Computational Settings

group water type 11 12
group sodium type 10

neighbor 2.0 bin
neigh_modify check yes

timestep 0.01

run_style verlet

velocity all create 300.0 713490 dist gaussian

#OPT 1:
fix 1 all temp/rescale 10 300.0 300.0 10.0 1.0
fix 2 water shake 0.0001 10 500 b 3 a 1

# End OPT 1

#OPT 2:
fix 3 all nve

fix 1 all temp/rescale 10 300.0 300.0 10.0 1.0
fix 2 water shake 0.0001 10 500 b 3 a 1
#OPT 2

#OPT 3:

fix 1 all nvt 300.0 300.0 0.001
fix 2 water shake 0.0001 10 500 b 3 a 1


# End OPT 3

thermo_style one
thermo 100

dump 1 all dcd 100 NaMMT_vsMD5ps.dcd
dump 2 all atom 100 NaMMT_vsMD5ps.atom
dump 3 all atom 100 NaMMT_vsMD5ps.xtc

dump_modify 2 scale no

restart 1000 restart_NaMMT_vsMD5ps.*

run 5000

dear aric,

thank you for your detailed description of the problems that you are seeing.

[...]

If I want to update the atom positions, it is advised to use a separate time
integration fix like fix nve (OPT 2). If I do include the time integration
fix nve, the atom positions do change and the temperature remains at the
desired level (300K). The SHAKE stat ave is constant and the delta changes
a bit from step to step. But the TotEng decreases (becomes more negative)
while the pressure increases (becomes more positive) and visually (via VMD)
the mineral collapses and LAMMPS reports and ERROR: Dump dcd of non-matching

here you have the indication of a serious problem!
it has nothing to do with the initial velocities.
more likely there is something not correct with
your potential parameters or you are not using
the correct potential types (pair/bond/angle/dihedral/improper styles).

# of atoms.
If I try to bypass the velocity-scaling and simply use fix nvt (OPT 3) to
try to set the temp, then the temperature, TotEng, and Pressure all increase
until the standard result is "nan" and eventually the run crashes. If I
comment out the velocity assignment while using OPT 3 and read in the last
restart file from OPT 1 (read_restart...), the result at step 1 is "nan" for
the SHAKE stats and for the themo output.

this is just a manifestation of the same issue. the nvt integrator
diverges, because your interactions are far too strong and then
you get too much kinetic energy in your system and then everything
goes haywire. it is worse with nvt compared to nve plus temp/rescale
because the latter can remove kinetic much more efficiently.

I would appreciate any advice, tips, or suggestions that could be provided
on how to assign the initial velocities for an MD run.

i would recommend you very carefully check your potential parameters
_and_ what would be the correct unit setting for lammps.

cheers,
     axel.

Dear Alex,

Thank you for your reply. You are correct, double-checking the potential parameters revealed that they were way off. My proliferation of test and input files resulted in me using the wrong read_data file. Also, and I don't know if it mattered, I failed to specify a mixing rule for the potential parameters (pair_modify mix...).

Thanks for taking the time to answer my question.

Aric