Moltemplate cannot create Data Atoms file

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

2 Likes