For this job we will use the cluster HYPATIA available at Empa. There we have access to parallel facilities with reserved nodes for the lecture. How to connect to HYPATIA:

Dear Student,

In order to be able to run simulations at high priority, today we will work on the Empa Cluster. We have created a personal account for you. Since the cluster is behind a firewall, we must connect to a gate machine (jumphost) to be allowed to access to the cluster.

Here the instructions to connect.

1) Decide a username for hypatia (it will be one between mmmstud01 and mmmstud25 ) 2) connect to the jumphost with the username (same for all) mmmstud

ssh -X mmmstud@jump1.empa.ch Password: will be communicated in the class

3) Connect to hypatia: ssh -X mmmstud02@hypatia (replace “02” by your number) password: same as username

4) you are in! 

[you@hypatia ~]$ mmm-init
[you@hypatia ~]$ cd /mnt/scratch/YOURUSER/
[you@hypatia ~]$ cp -r /mnt/scratch/psd/exercise_12 .
[you@hypatia ~]$ cd exercise_12

Reproducing a PMF calculation for adsorption of an organic molecule on KCl

We will today refer to the paper available at 10.1021/acs.jpcc.5b12028 about the entropic effects of adsorption of pretty large molecules on KCl. We will use the software LAMMPS for performing a series of MD simulations at fixed molecular height above the surface. Then we will use integration of the average force between the molecule and the substrate, along the adsorption path, to extract the free energy difference between two configurations: the “very far” configuration and a configuration along the path.

  • After you copied the exercise_12 directory and entered it, look at the mol_sub.xyz file (by editing or vmd) and understand the geometry of the system: which range of distances should you consider for MD runs?
  • Create a directory T_300 and copy into it the following files and cd into the directory:
 cp run* *pot c1.topo md_temp.inp get_pot_mean_force T_300  ; cd T_300
  • At this point edit the run_distance_loop script and insert the list of distances you want to simulate.
  • How many steps will you run per each distance? This will be decided by editing the input file.
md_temp.inp 
# PMF Input script template
# (1) Initialisation
units real
atom_style full
boundary p p p
pair_style hybrid/overlay lj/charmm/coul/long 10.0 12.0 buck/coul/long 12.0 morse 10.0
bond_style harmonic
angle_style harmonic
dihedral_style charmm
kspace_style pppm 0.0001
# (2) Define atoms
read_data c1.topo
group molecule molecule 1
group substrate molecule 2
group bottom molecule 3
# (3) Settings
# freeze bottom layer of substrate to prevent it from drifting
fix 2 bottom setforce 0 0 0
include new_defpot.pot
include best_all.pot
neighbor 2.0 bin
neigh_modify delay 0
timestep 1.0
# (4) NVT Dynamics
fix temp1 molecule nvt temp _TEMP_ _TEMP_ 100
fix temp2 substrate nvt temp _TEMP_ _TEMP_ 100
velocity                        all create _TEMP_ 293288
fix PMF molecule recenter NULL _Y_ NULL
compute temp_molecule molecule temp
compute yforce molecule group/group substrate kspace yes
thermo_style custom time c_yforce[2] etotal pe c_temp_molecule temp ke evdwl press enthalpy
thermo_modify flush yes
thermo 50
dump xyz all xyz 100000000 mol_sub.xyz
dump_modify xyz element C C H C N K Cl
dump coord all dcd 5000 trajectory.dcd
restart _NSTEPS_ TCB_PMF.restart
run _NSTEPS_
  1. Edit the input file to run 50 picoseconds, to write the restart at the end, to thermalize at the wished temperature
  2. REPLACE _NSTEPS_ WITH THE NUMBER OF STEPS AND _TEMP_ WITH THE DESIRED TEMPERATURE IN THE md_temp.inp FILE
  3. The length of the simulation is small to get converged averages, but the run will thus last a few minutes per distance making the exercise feasible on 16 cores.
  4. Now run the chain of simulations
 qsub run_distance_loop
  1. Directories R_<distance> will be created (check with ls -ltr )
  2. You can enter a certain directory (e.g., R_12.5), check the log.lammps for the evolution of a trajectory, and visualize the trajectory with: vmd mol_sub.xyz trajectory.dcd
  3. At the end of everything, you can perform averages and integrals by understanding (and using) the script get_pot_mean_force FROM WITHIN THE DIRECTORY T_300:
./get_pot_mean_force
  1. THIS PRODUCES THE POTENTIAL OF MEAN FORCE pot_mean_force for that temperature
  2. Repeat the same for other temperatures (e.g.) 10 K and 800 K.
  3. Potentials can be visualized by adapting the gnuplot script in the exercise_12 directory: pot_mean_force.gnu. Edit and adapt, then:
 gnuplot
gnuplot> load "pot_mean_force.gnu"
Check the already finished runs in REF_300: 
cd REF_300
cp ../get_pot_mean_force .
./get_pot_mean_force
gnuplot
gnuplot > plot "pot_mean_force" w lp
  1. Report adsorption (free) energy and equilibrium distance for the three temperatures. What do you note?
  2. Compare the results with the paper
  3. Compare your results with results extracted from REF_10, REF_300, REF_800 that contain longer equilibrations and averaging (about 0.5-1 ns each distance). Where do you get more differences? Why?
  4. Can you say something about the entropy contributions?