NPT molecular dynamics in ASE

Dear ASE developers and community,

I am using ASE to conduct MD simulation of a polymer system, using ANI calculator. I am loading a structure from a LAMMPS data file using:

atoms = read(‘Poly.data’, format=‘lammps-data’, Z_of_type={1: 6, 2: 17, 3: 1, 4: 7})

T = 300
MaxwellBoltzmannDistribution(atoms, temperature_K=T)

calc_ani2x = torchani.models.ANI2x().ase() # import ANI-2x
atoms.set_calculator(calc_ani2x) # Specify ANI-2x as the calculator

print(“Beginning NVT equilibration…”)

dyn1 = NVTBerendsen(atoms, trajectory= ‘md.traj’, timestep=1 * units.fs, temperature_K=300, taut=0.5 * 100 * units.fs)

logger_nvt = MDLogger(dyn1, atoms, ‘md.log’, header=True, stress=False, mode=“w”)
dyn1.attach(logger_nvt, interval=10)
dyn1.run(100)

print(“Beginning NPT equilibration…”)

pressure_in_bar = 1.0
pressure_in_ev_per_ang3 = pressure_in_bar / 160217.66208

dyn2 = NPT(atoms, timestep=1 * units.fs, temperature_K=300, ttime=25 * units.fs, externalstress=pressure_in_ev_per_ang3, pfactor= 30 * units.fs, mask=(1, 1, 1),
trajectory= ‘md.traj’, append_trajectory=True)

logger_npt = MDLogger(dyn2, atoms, ‘md.log’, header=True, stress=False, mode=“a”)
dyn2.attach(logger_npt, interval=10)
dyn2.run(100)

However, when I check the process of the md simulation, the volume of the simulation box does not change much, could somebody provide some information on how to choose values of ttime and pfactor to yield meaningful simulation results? My system size is not that big, just hundreds of atoms.

Thank you all. Appreciate it!

1 Like

Hello, any update so far?
Have you tried minimizing the energy before starting the equilibration? What about boundary conditions? Are they periodical?