Generation of Water in Water Trajectories

First, we need an input coordinate file in TINKER XYZ (or PDB) format

                    648 Water Cubic Box (18.643 Ang, 216 AMOEBA)
                    1 O 8.039430 5.868035 0.492777 1 2 3
                    2 H 7.581494 5.022666 0.403759 2 1
                    3 H 8.287056 6.062440 -0.397074 2 1
                    4 O 0.114840 -8.876599 6.445852 1 5 6
                    5 H 0.958925 -8.446494 6.405247 2 4
                    6 H 0.247235 -9.794471 6.137893 2 4
                    7 O -6.576367 -0.252600 8.104632 1 8 9
                    8 H -6.645499 0.689925 7.896356 2 7
                    9 H -6.684537 -0.400760 9.103747 2 7
                

Next, is the property file (like a TINKER keyword file). Note that Force Field X stores parameter files internally.

                    forcefield amoeba-water
                    a-axis 18.643
                    spacegroup P1
                    polar-eps 0.01
                

Finally, we use a Force Field X script to run the sampling.

                    
// Coordinate file to open (can be XYZ or PDB)
String fileName = "examples/watersmall.xyz";

// Beginning of the ligand atom range.
int ligandStart = 1;

// End of the ligand atom range.
int ligandStop = 3;

// Number of electrostatics lambda windows.
int elecWindows = 10;

// Number of soft core van der Waals lambda windows.
int vdwWindows = 10;

// Number of equilibration MD steps.
int eSteps = 100000;

// Number of MD steps per window.
int nSteps = 100000;

// Time step in femtoseconds.
double timeStep = 1.0;

// Frequency to print out thermodynamics information in picoseconds.
double printInterval = 0.01;

// Frequency to save out coordinates in picoseconds.
double saveInterval = 0.1;

// Temperature in degrees Kelvin.
double temperature = 300.0;

// Things below this line normally do not need to be modified.
// =============================================================================

import org.apache.commons.io.FilenameUtils;

import ffx.algorithms.MolecularDynamics;
import ffx.algorithms.dynamics.thermostats.Thermostat.Thermostats;

import ffx.potential.PotentialEnergy;
import ffx.potential.bonded.Atom;

// Open the system
open(fileName);
String lambdaName = FilenameUtils.removeExtension(fileName);

// Apply the ligand atom selection
Atom[] atoms = active.getAtomArray();
for (int i = ligandStart; i <= ligandStop; i++) {
    Atom ai = atoms[i - 1];
    ai.setApplyLambda(true);
}

// Select an availble thermostat [ BUSSI, BERENDSEN, ISOTHERMAL ]
Thermostats thermostat = Thermostats.BUSSI;

// Create a MolecularDynamics instance
MolecularDynamics molDyn = new MolecularDynamics(active, null, thermostat);

// Equilibrate the system
boolean initVelocities = true;
molDyn.dynamic(eSteps, timeStep, printInterval, -1, temperature, initVelocities);

// Retrieve the PotentialEnergy instance
PotentialEnergy energy = active.getPotentialEnergy();

// Turn off the polarizability and multipoles for the selected atoms
for (int i=0; i <= elecWindows; i++) {
    double lambda = 1.0 - (double) i / (double) elecWindows;
    energy.setElectrostaticsLambda(lambda);
    molDyn.setArchiveFile(new File(lambdaName + "_elec_" + String.format("%5.3f", lambda) + ".arc"));
    molDyn.dynamic(nSteps, timeStep, printInterval, saveInterval, temperature, initVelocities);
}

// Turn off the vdW potential for the selected atoms
for (int i=1; i <= vdwWindows; i++) {
    double lambda = 1.0 - (double) i / (double) vdwWindows;
    energy.setSoftCoreLambda(lambda);
    molDyn.setArchiveFile(new File(lambdaName + "_vdw_" + String.format("%5.3f", lambda) + ".arc"));
    molDyn.dynamic(nSteps, timeStep, printInterval, saveInterval, temperature, initVelocities);
}