# Import the force field.
import /usr/local/moltemplate/moltemplate/force_fields/oplsaa.lt
import formamide.lt # after oplsaa.lt, as it depends on it.

# Create the input sample.
solv = new _FAM [5].move( 4.6, 0, 0)
                [5].move( 0, 4.6, 0)
                [5].move( 0, 0, 4.6)
solv[*][*][*].move(-11.5, -11.5, -11.5)

# Set the simulation box.
write_once("Data Boundary") {
 -11.5 11.5 xlo xhi
 -11.5 11.5 ylo yhi
 -11.5 11.5 zlo zhi
}

# Create an input deck for LAMMPS.
write_once("In Init"){
 # Input variables.
 variable run    string solv_01   # output name
 variable ts     equal  1         # timestep
 variable temp   equal  300       # equilibrium temperature
 variable p      equal  1.        # equilibrium pressure
 variable d      equal  1000      # output frequency
 variable equi   equal  5000      # Equilibration steps
 variable prod   equal  30000     # Production steps

 # PBC (set them before the creation of the box).
 boundary p p p
}

# 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\} iso \$p \$p \$\{pcouple\}
 velocity all create \$\{temp\} 858096 dist gaussian
 # Short runs to update the PPPM settings as the box shinks.
 run    \$\{equi\} post no
 run    \$\{equi\} post no
 run    \$\{equi\} post no
 run    \$\{equi\}
 # From now on, the density shouldn't change too much.
 run    \$\{prod\}
 unfix NPT
}