Disclaimer: this thread is slightly off-topic in this forum, but it’s still helpful to teach how to use Moltemplate to create suitable initial geometries for MD simulations and a complete input deck for LAMMPS.
Fabrizio, your idea was correct. However, you made many syntactic mistakes:
- The files should be plain text, but yours seem to be formatted with a word editor: quotes and hyphens are replaced with special characters, e.g.
wrong:write_once(“In Settings”) {
correct:write_once("In Settings") {
or
wrong:angle_coeff @\angle:HC–CT–CZ 35.0 108.5
correct:angle_coeff @angle:HC-CT-CZ 35.0 108.5
- Moltemplate variables should be declared without escaping characters:
wrong:@\atom:CT 12.010700
correct:@atom:CT 12.010700
- In the molecule definition (file
acn.lt
), the bond list must be defined with the variable used in the"Data Atoms"
section:
wrong:$bond:ch1 $atom:CT $atom:HC_1
correct:$bond:ch1 $atom:c2 $atom:h21
- The geometry that you specified has a lot of overlaps. But without creating the data file, it’s hard to tell!
I rearranged your input files more logically. This should get you started. The force field does not contain any dihedral or improper term to control the rotation of the terminal CH3 group, you may want to add one.
Here are the revised files.
File: ForceField_ACN.lt
ForceField {
write_once("Data Masses") {
@atom:CT 12.010700
@atom:CZ 12.010700
@atom:NZ 14.006700
@atom:HC 1.007940
}
write_once("In Settings") {
pair_coeff @atom:CT @atom:CT 0.066 3.30
pair_coeff @atom:HC @atom:HC 0.015 2.50
pair_coeff @atom:NZ @atom:NZ 0.170 3.20
pair_coeff @atom:CZ @atom:CZ 0.066 3.30
}
# Bonds
write_once("Data Bonds By Type") {
@bond:CT-CZ @atom:CT @atom:CZ
@bond:CT-HC @atom:CT @atom:HC
@bond:CZ-NZ @atom:CZ @atom:NZ
}
write_once("In Settings") {
bond_coeff @bond:CT-CZ 385.0 1.458
bond_coeff @bond:CT-HC 410.0 1.087
bond_coeff @bond:CZ-NZ 650.0 1.157
}
# Angles
write_once("Data Angles By Type") {
@angle:HC-CT-CZ @atom:HC @atom:CT @atom:CZ
@angle:CT-CZ-NZ @atom:CT @atom:CZ @atom:NZ
@angle:HC-CT-HC @atom:HC @atom:CT @atom:HC
}
write_once("In Settings") {
angle_coeff @angle:HC-CT-CZ 35.0 108.5
angle_coeff @angle:CT-CZ-NZ 150.0 180.0
angle_coeff @angle:HC-CT-HC 35.0 109.5
}
# You may want to add dihedral and improper angles too!
# These settings are needed for every simulation, so it's good
# to store them along with the force field.
write_once("In Init") {
units real
atom_style full
pair_style lj/cut/coul/long 11.0 11.0
bond_style harmonic
angle_style harmonic
dihedral_style opls
improper_style harmonic
pair_modify mix geometric
special_bonds lj/coul 0.0 0.0 0.5
kspace_style pppm 0.0001
}
}
File: acn.lt
import ForceField_ACN.lt
ACN inherits ForceField {
write("Data Atoms") {
$atom:c2 $mol:acn @atom:CT 0.080 -0.527520 -0.736650 -0.061930
$atom:c1 $mol:acn @atom:CZ 0.460 -0.490180 0.691350 0.123090
$atom:n1 $mol:acn @atom:NZ -0.560 -0.460210 1.840750 0.272090
$atom:h21 $mol:acn @atom:HC 0.060 0.490830 -1.139500 -0.140720
$atom:h22 $mol:acn @atom:HC 0.060 -1.074560 -0.989770 -0.979630
$atom:h23 $mol:acn @atom:HC 0.060 -1.028760 -1.219740 0.787100
}
# A list of atoms must be defined using the variable names used in the previous section.
write("Data Bond List") {
$bond:cc $atom:c1 $atom:c2
$bond:ch1 $atom:c2 $atom:h21
$bond:ch2 $atom:c2 $atom:h22
$bond:ch3 $atom:c2 $atom:h23
$bond:cn $atom:c1 $atom:n1
}
# The group command can go in the same file as the rattle command.
write_once("In Constraints") {
group acn type @atom:CT @atom:CZ @atom:NZ @atom:HC
fix fRattleAcn acn rattle 0.0001 10 100 b @bond:CT-HC a @angle:HC-CT-HC
}
}
File: acn_box.lt
# I use moltemplate in a slightly different way,
# with a master file containing all the commands needed
# to create a valid input deck for LAMMPS.
# Define the input variables.
write_once("In Init"){
# Input variables.
variable run string acn_box # output log
variable ts equal 1 # timestep
variable temp equal 300 # equilibrium temperature
variable p equal 1. # equilibrium pressure
variable d equal 1000 # output frequency
variable prod equal 30000 # Production steps
# PBC (set them before the creation of the box).
boundary p p p
}
# Import the force field.
import acn.lt
# Create the input sample.
# Changing geometry to avoid bad contacts.
solv = new ACN [7].move( 5, 0, 0)
[7].move( 0, 5, 0)
[4].move( 0, 0, 5)
solv[*][*][*].move(-17.5, -17.5, -10)
# Set the simulation box.
write_once("Data Boundary") {
-17.5 17.5 xlo xhi
-17.5 17.5 ylo yhi
-10 10 zlo zhi
}
# Run an NPT simulation.
write_once("In Run"){
# Derived variables.
variable tcouple equal \$\{ts\}*100
variable pcouple equal \$\{ts\}*1000
# Output.
thermo \$d
thermo_style custom step etotal evdwl ecoul elong ebond eangle edihed eimp &
ke pe temp press vol density cpu
thermo_modify flush yes
# Trajectory.
dump TRJ all dcd \$d \$\{run\}.dcd
dump_modify TRJ unwrap yes
# Thermalisation and relaxation, NPT ensemble.
timestep \$\{ts\}
fix NPT all npt temp \$\{temp\} \$\{temp\} \$\{tcouple\} aniso \$p \$p \$\{pcouple\}
velocity all create \$\{temp\} 858096 dist gaussian
run \$\{prod\}
unfix NPT
}
Note that here, I have used the escape character to protect LAMMPS variables inside the Moltemplate script, e.g. \$\{run\}
.
Compile with:
moltemplate.sh acn_box.lt -overlay-all
Run with:
mpirun -np 4 lmp_2Aug23 -in acn_box.in -l acn_box.log