Free energy perturbation module

Solvation energy of group IIA ions

Objective: to create a model of a solvated divalent ion and study the change in solvation energy by the free energy perturbation technique.

Reference: Bartolotti, L.J., Pedersen, L. G., Charifson, P. S. "Long-range nonbonded attractive constants for some charged atoms" Journal of Computational Chemistry (1991), 12, 1125-1128

Criteria for analysis:

    1. Relative free energy of solvation of different ions in the series Mg2+, Ca2+, Sr2+.
    2. Radial distribution function of the solvent surrounding the ion.

Option 1: Setup your own car and mdf files.

Option 2: Download the car and mdf files.

ca_to_sr.car

ca_to_sr.mdf

sr_to_ca.car

sr_to_ca.mdf

mg_to_ca.car

mg_to_ca.mdf

ca_to_mg.car

ca_to_mg.mdf

Before we proceed with the energy minimization and dynamics we must select a forcefield. The default forcefield in DISCOVER is CVFF (the consistent valence forcefield which is also default forcefield in DISCOVER3). This forcefield does not have parameters for Mg, Ca, or Sr. These parameters have been derived by Bartolotti [J. Comp. Chem., (1991), 12, 1125-1128]. The parameters have been incorporated into a forcefield file cvff_ion.frc given below. You will need to copy this forcefield and the associated cvff_ion.rlb and cvff_ion_templates.dat files to the directory where you are going to run DISCOVER.

cvff_ion.frc

cvff_ion.rlb

cvff_ion_templates.dat

To compile the force field and create a binary cvff_ion.bin file you can use the script run_binff

You can make the script executable by typing: chmod +x run_binff at the UNIX prompt

If you are following option 1 you may need to use the forcefield definitions to modify the structure you are creating for your car and mdf files.

You can obtain input files in two ways.

Option 1: generate your own input files.

Option 2: download the input files from the list below.

sr_to_ca.inp

ca_to_sr.inp

mg_to_ca.inp

The stages to be run are:

  1. energy minimization
  2. 10 steps of equilibration and dynamics

The steps of equilibration and dynamics are exactly the same as those you ran for either the argon or water simulations, however the set up will be different. In the free energy perturbation calculation the equilibration and dynamics are carried out in 10 or more steps as the parameters are changed from those of ion 1 to ion 2.

FOR THIS MODULE YOU WILL WANT TO USE THE VERSION KNOWN AS DISCOVER (DISCOVER2.9) AND NOT THE DISCOVER_3 (DISCOVER3.0) PROGRAM USED FOR THE PREVIOUS MODULES.

To run discover interactively, follow the procedure below.

The first prompt will be

>> ENTER to file name of PREFIX (<RETURN>=QUIT)

>> [Enter the root name here. Be sure that the inp, car, and mdf files for minimization all have the same root name]

Next you are prompted for the forcefield

>> ENTER 1, 2, 3, 4 or 5, (<RETURN> = 1) 5

>> ENTER the potential library name (<RETURN> = QUIT) cvff_ion.bin

>> ENTER the nice number for running discover as a background process.

The default is 19.

>> [Hit return to use the default]

NOTE: IN THE UNIX OPERATING SYSTEM PRIORITY FOR JOB GOES FROM 1 HIGHEST TO SOME LOW NUMBER. BY GIVING YOURSELF A NICE LOW NUMBER YOU ALLOWING YOUR JOB TO RUN WITHOUT INTERRUPTING HIGHER PRIORITY INTERACTIVE JOBS LIKE YOUR OWN INSIGHT SESSION!

>>Enter the number of processors – default is 32

>> 1

NOTE: THIS PROGRAM SEEMS TO RUN SLOWER IN MULTIPROCESSOR MODE. THIS PROBABLY HAS TO DO WITH ALLOCATION OF THE PROCESSORS ON THE ORIGIN 2000. DO NOT ATTEMPT TO RUN WITH MULTIPLE PROCESSORS.

>> Do you wish to start discover now?

The default is YES

>> <Return>

Now discover will run for about 10 minutes for minimization and 25 minutes for the free energy perturbation.

Before running the dynamics stage for free energy perturbation be sure to rename the cor file from minimization (final coordinates of the energy minimized solvent box) to be the car file for the subsequent free energy perturbation calculation.

To run the free energy perturbation calculation you need to use the command discover. Here the input will be interactive. It will have exactly the sequence of steps listed above for the minimization.

To examine the progress of the job observe the *.pek file

Analysis

Once the job is complete you can analysis the results using the DECIPHER module.

Read in the cor file using the Molecule/Get command in INSIGHT. After running dynamics you will see an apparently unstructured argon solvent as shown below.

Use the side pull-down menu to run DECIPHER. In decipher use the Configurations/Get command to generate a menu. Use the Files list on the right-hand side to get the his file.

Configurations File [history file]

Assembly/molecule [ASSEMBLY_NAME]

Load_Cell (check the box to turn it yellow)

Load_Pressure (do not check the box)

Load_Velocities (check the box to turn it yellow)

Frame Spec [*] (reads all of the frames)

Select execute.

To calculate the radial distribution function averaged over the entire trajectory (i.e. averaging all of the frames in the history file) use the

Spectra/RadialDistFunc command.

Distribution_type [pair]

Molecule_Spec [OHHC_B0]

Molecule_2_Spec [OHHC_B0_CELL_SOLVEA]

Coord_Number (check the box to turn it yellow)

Bulk_Density (check the box to turn it yellow)

Data_Sources [trajectory] (check the box to turn it yellow)

Background (check the box to turn it yellow)

Create_TBL_File [rdf] (check the box to turn it yellow)

To calculate the mean squared displacement use the Spectra/Mean_Sq_Disp command.

Molecule_Spec [OHHC_B0]

Data_Source [trajectory]

Frame_Spec [*]

Output

MSD (check the box to turn it yellow)

Dist_Traveled (check the box to turn it yellow)

Background (check the box to turn it yellow)

Create_TBL_File [msd] (check the box to turn it yellow)