Build and equilibrate the system

The system consists of a single Na+ ion in water.  To build it, start with a .gro file with a single Na+ ion located at the center of a 2.5×2.5×2.5nm box.  Use gmx solvate to add water.  Equilibrate with a standard protocol of energy minimization, followed by 100 ps each of NVT and then NPT simulation, with a timestep of 2fs.  A single core is plenty for such a small system.    Make the production run 1ns; set nsfxout and nsftout to write 1000 frames.  The position of the Na+ ion should be frozen throughout, with

freezegrps = Ion
freezedim = Y Y Y

Add a test charge

To add a test charge, we alter the .top file to change the definition of a selected water molecule.  Copy the file oplsaa.ff/spce.itp as spceq.itp; change the molecule name to SOLq; add +1 charge to the oxygen.  Finally, remove the SETTLES declarations (only one set of SETTLES are permitted, and these are already used in spce.itp). 

Make a .top file systemq0.top, which #includes “spceq.itp” after spce.itp.  The [ molecules ] declaration should read:

[ molecules ]
; Compound         #mols
NA                         1
SOL                       NTOP
SOLq                     1
SOL                       NBOT

This is a “template” .top file; we will replace the values NTOP and NBOT from within the job script, using the Unix utility sed (Stream EDitor), described below.  So we can access the test molecule as a group, add the following “template lines” to the bottom of an index file systemq.ndx:

[ testq ]
OW HW1 HW2

OW, HW1, and HW2 will likewise be replaced using sed in the job script with the atom numbers of the selected water molecule. 

Loop over molecules

The job script “testAllq.sh” consists of a loop that does the following tasks:

•  select a water molecule at random
• use sed to create .top file from template systemq.top
• use sed to create .ndx file from template systemq.ndx
• rerun existing trajectory without test charge; write force on selected molecule to file
• rerun trajectory with test charge; write selected force, positions to file
• remove backup files on each iteration with rm “#”*”#”

Several tricks and techniques are employed in the brief resulting script (my version is 23 lines of code).  Here are some snippets to serve as a guide (the full script is included with all necessary input files in the tarball “screening.tar”):

pass variable to script

Use a variable KMAX to limit the loop.  Pass its value to the script when the job is submitted, with “qsub -v KMAX=500…”

“while” loop

A “while” loop in bash scripts has the form

let K=1
while (($K < $KMAX))
do
…
K=$((K+1))
done

random numbers

Each time the system variable RANDOM is accessed, it returns a random integer from 0 to 32767 inclusive.  To generate a random number R between X and Y inclusive, use

BASE=$((Y-X+1))
R=$(($(($RANDOM%$BASE))+X))

(Here “%” is the “modulo” operator, i.e., $A%$B gives the remainder on dividing A by B.)

If you have NTOTAL water molecules, you want NTOP to be between 1 and NTOTAL-2, and NBOT to be (NTOTAL-1-NTOP) [this way, both NTOP and NBOT are 1 or larger].  Once you have NTOP, the selected molecule is number NTOP+1; its atoms OW, HW1, HW2 are numbers 1+3*NTOP+1, 1+3*NTOP+2, 1+3*NTOP+3 (the first “1” is the Na+ ion, the next “3*NTOP” are the atoms of previous water molecules).

templates with sed

An example of how to use sed to fill in a “template file” is as follows:

sed “s/NTOP/${NTOP}/;s/NBOT/${NBOT}/” < $SYSTEM/systemq0.top > $SYSTEM/systemq.top

Here $SYSTEM is the system folder (containing .top, .itp, .ndx files); NTOP and NBOT are the previously computed numbers of waters above and below the selected molecule.

sed does substitutions on each line it reads from stdin (STanDard INput), and writes to stdout (STanDard OUTput), here substituting the script variable ${NTOP} for the literal string “NTOP”, and ${NBOT} for NBOT.  Both substitutions are done (in order) on each line sed processes.

-rerun

An example of how to rerun an existing trajectory is as follows:

gmx grompp -f $MDP/run.mdp -c $INIT -p $SYSTEM/systemq.top -n systemq.ndx -o rerunq.tpr -maxwarn 20
gmx mdrun -nt 1 -deffnm rerunq -rerun run.trr

Here $INIT is the initial frame used in the original run, which produced a trajectory run.trr. systemq.top is the .top file filled in by sed from the template, and systemq.ndx likewise the filled-in .ndx file.  The flag -rerun tells mdrun to rerun the trajectory run.trr:  that is, to compute forces and energies using the new .tpr file, but move the atoms precisely as in the existing trajectory run.trr. 

write forces, positions

The Gromacs utility gmx traj can write positions or forces of selected atoms to an .xvg file.  Here are examples of how to use it:

echo “7” | gmx traj -f rerunq.trr -s rerunq.tpr -n systemq.ndx -com -of forceq_${K}.xvg

The construction “echo “7” |” effectively “types” input (here selecting the group to be output) which is then “piped” with | to gmx traj, thus allowing us to “script” an interactive program.  The output file name is constructed from the loop variable.  Flag -com writes the force for the COM of the group; -of indicates output of forces.  To output positions of all atoms in a group (as for the selected water molecule), we use 

echo “7” | gmx traj -f rerunq.trr -s rerunq.tpr -n systemq.ndx -ox posn_${K}.xvg

[Page 3:  analyze the results.]