User Tools

Site Tools


exercises:2015_ethz_mmm:c2h2_bond_energy

This is an old revision of the document!


C2H2 and C2H4 bond energy

TO USE THE FUNCTION LIBRARY (VERSION UP TO DATE IN DANIELE PASSERONE DIR) IN THE INTERACTIVE SHELL: you@eulerX ~$ . ~danielep/Scripts/m_functions.bash
<note warning>
Before you start it is **ABSOLUTELY NEСESSARY** to update the function library. So from you $HOME directory please run:
<code>
you@eulerX ~$ ./update_functionlibrary

Pay attention to the output of this program. It may give you some helpful information

If you don't have the file “update_functionlibrary” please download it from here update_functionlibrary.zip and unzip into your $HOME directory

you@eulerX ~$ unzip update_functionlibrary.zip
you@eulerX ~$ chmod 755 update_functionlibrary

Or if you are working on the remote computer you can also do the following:

you@eulerX ~$ wget http://www.cp2k.org/_media/exercises:2015_ethz_mmm:update_functionlibrary.zip
you@eulerX ~$ unzip exercises:2015_ethz_mmm:update_functionlibrary.zip
you@eulerX ~$ chmod 755 update_functionlibrary

</note>

REMEMBER: this is the command to load the module for the cp2k program:
you@eulerX ~$ module load new cp2k

and to submit the job:

you@eulerX ~$ bsub < jobname

Download the 2.1 exercise into your $HOME folder and unzip it.

you@eulerX ~$ wget http://www.cp2k.org/_media/exercises:2015_ethz_mmm:exercise_2.1.zip
you@eulerX ~$ unzip exercises:2015_ethz_mmm:exercise_2.1.zip
All files of this exercise (input and scripts are all commented) can be downloaded from the wiki: exercise_2.1.zip

In this exercise you will test the general AMBER force field, which is suitable for a large class of molecules. In this example you will compute the classical energy of a simple acetylene molecule and apply small variations to the C-C bond lengths.

Then go to the directory “exercise_2.1/c2h2/”.

you@eulerX ~$ cd exercise_2.1/c2h2/

The input file structure of the template is the following:

c2h2.inp.templ
&FORCE_EVAL                     ! This section defines method for calculating energy and forces
  METHOD FIST                   ! Using Molecular Mechanics
  &MM
    &FORCEFIELD                 ! This section specifies forcefield parameters
      parm_file_name c2h2.top
! This file contains force field parameters
      parmtype AMBER
      &SPLINE                   ! This section specifies parameters to set up the splines used in the nonboned interactions (both pair body potential and many body potential)
        EMAX_SPLINE 15.0        ! Specify the maximum value of the potential up to which splines will be constructed
      &END
    &END FORCEFIELD
     &POISSON                   ! This section specifies parameters for the Poisson solver
      &EWALD                    ! This section specifies parameters for the EWALD summation method (for the electrostatics)
        EWALD_TYPE SPME         ! Smooth particle mesh using beta-Euler splines
        ALPHA .36
        GMAX 128                ! Number of grid points
      &END EWALD
    &END POISSON
    &PRINT                      ! This section controls printing options
      &FF_INFO
      &END
      &FF_PARAMETER_FILE
      &END
    &END
  &END MM
  &SUBSYS                       ! This section defines the system
    &CELL                       ! Unit cell set up
      ABC 50 50 50              ! Lengths of the cell vectors A, B, and C
    &END CELL
    &COORD                      ! This section specify all the atoms and their coordinates
H                   0 0 0 
C                   _C1_  0 0 
C                   _C2_  0 0
H                   _H2_  0 0 
     &END COORD
    &TOPOLOGY
      CHARGE_BETA               ! Read MM charges from the BETA field of PDB file
      CONNECTIVITY AMBER        ! Use AMBER topology file for reading connectivity
      CONN_FILE_NAME c2h2.top
! This file contains connectivity data
    &END TOPOLOGY
    &PRINT
      &TOPOLOGY_INFO
        AMBER_INFO
      &END
    &END
  &END SUBSYS
&END FORCE_EVAL

&GLOBAL                         ! Section with general information regarding which kind of simulation to perform an parameters for the whole PROGRAM
  PRINT_LEVEL LOW               ! Global print level
  PROJECT c2h2                  ! Name of the project. This word will appear as part of a name of all ouput files (except main ouput file, specified with -o option)
  RUN_TYPE ENERGY               ! Calculating of energy for the fixed position of atoms
&END GLOBAL                                                                                                                                                                                

The script to modify the C-C distance in the input and generate the curve is the following:

c2h2.chain
#
# task 2.1 C2H2 bond strength
#
# the distances of equilibrium
# and increment delta
#
dist_cc0=1.181
dist_hc=1.066
delta=0.0002
 
#
# loads the functions
# this script uses the functions: change replace getcolumn add list
#
module load new cp2k
. ~/Scripts/myfunctions.bash
 
 
#
# deletes old output file
#
rm curve.amber.c2h2
 
#
# loop over all values of n
#
 
for n in `list -5 5`
do
  dist_cc=`change $dist_cc0 $n $delta`
  c2=`add $dist_hc $dist_cc `
  h2=`add $c2 $dist_hc `
  #
  # replaces coordinates in the input
  #
  replace _C1_ $dist_hc < c2h2.amber.inp.templ | replace _C2_ $c2  | replace _H2_ $h2 > c2h2.$dist_cc.inp
  #
  # Runs cp2k
  #
  echo RUNNING cp2k with the input c2h2.$dist_cc.inp
  # bsub -n 1 mpirun cp2k.popt -i c2h2.$dist_cc.inp > c2h2.$dist_cc.out
  cp2k.popt -i c2h2.$dist_cc.inp > c2h2.$dist_cc.out
  #
  # gets the energy from the output (in a.u.)
  #
  energy=` getcolumn 'ENERGY|' 9 < c2h2.$dist_cc.out `
  echo $dist_cc $energy  >> curve.amber.c2h2
#
# loop over n
#
done
 
#
# Cleaning...
#
mkdir Logs
mv c2h2.*.inp c2h2.*.out Logs
rm RUN* *restart*

At this point submit the job grid, first loading the module for cp2k entering

you@eulerX c2h2$ module load new cp2k
you@eulerX c2h2$ bsub < c2h2.chain

Finally, the gnuplot code to fit the curve to a parabolic function - and to extract the parameters of the bond interaction is the following

you@eulerX c2h2$ gnuplot fit.gnu
fit.gnu
#!/usr/bin/gnuplot
f(x)=k*(x-x0)*(x-x0)+e0         # definition of the function to be fitted
delta = 0.002                   #
e0=0.007720330647107            # initial value of e0
k = 500                         # initial value of k
x0=1                            # initial value of x0
#a0 = 1.3327436840              # 
set xlabel 'd (Angstrom)'       # label of the x-axis
set ylabel 'E (kcal/mol)'       # label of the y-axis
set title "c2h2 triple bond"    # main title

fit f(x) 'curve.a.c2h2' u ($1):(($2)*627) via k,x0,e0
# fitting procedure
#    ^ function to be fitted
#             ^ take data from the file curve.a.c2h2
#               and use the first column as X and the second one
#               ax Y (all values of the second column are multiplied by 627. 1 Hartree = 627 kcal/mol)
#                                               ^ k, x0, e0 - are parameters to be adjusted

plot 'curve.amber.c2h2' u ($1):(($2)*627)-e0 t 'computed AMBER' w lp lw 2 ps 2 pt 7, f(x)-e0 t 'k*(x-x0)^2' w l lw 2
#         ^ take data for the first plot from the file curve.a.c2h2
#                               ^ use the firts column as X and the second as Y
#                                               ^ title of the first plot: computed AMBER
#                                                              ^ draw points and lines
#                                                                ^ line width = 2 
#                                                                      ^ point size = 2
#                                                                           ^ point type = 7 
#                                                                                  ^ plot function f(x)-e0
#                                                                                               ^ title of the second plot: k*(x-x0)^2
#                                                                                                        ^ draw only line
#                                                                                                             ^ line width = 2


pause mouse                     # exit by right-clicking on the plot window
#    EOF

Compare the values that you obtain with the ones listed in the “human readable” potential file c2h2-force_field.pot that was generated by cp2k.

Now, perform the same exercise in another directory for the molecule C2H4.

you@eulerX c2h2$ cd ../c2h4
  • The new input c2h4.inp.templ is similar to the previous one, only with different project name (C2H4) and coordinates:
H   0 3.06 0 
H   0 4.94 0 
C   _C1_ 4.0 0
C   _C2_ 4.0 0
H   _H3_ 3.06 0
H   _H3_ 4.94 0 
  • The script to modify the bond lengths changes slightly as well
c2h4.chain
#
# task 2.2 C2H4 bond strength
#
# the initial deltax of  C, C, H, H
# and increment delta
#
dist_hcx=0.560
dist_cc0=1.324
delta=0.002
 
#
# loads the functions
# this script uses the functions: change replace getcolumn add list
#
module load new cp2k
. ~/Scripts/myfunctions.bash
 
 
#
# deletes old output file
#
rm curve.amber.c2h4
 
#
# loop over all values of n
#
 
for n in `list -5 5`
do
  dist_cc=`change $dist_cc0 $n $delta`
  c2=`add $dist_hcx $dist_cc`
  h3=`add $c2 $dist_hcx `
  #
  # replaces coordinates in the input
  #
  replace _C1_ $dist_hcx < c2h4.amber.inp.templ | replace _C2_ $c2  | replace _H3_ $h3 > c2h4.$dist_cc.inp
  #
  # Runs cp2k
  #
  echo RUNNING cp2k with the input c2h4.$dist_cc.inp
  cp2k.popt -i c2h4.$dist_cc.inp > c2h4.$dist_cc.out
  #
  # gets the energy from the output (in a.u.)
  #
  energy=` getcolumn 'ENERGY|' 9 < c2h4.$dist_cc.out ` 
  echo $dist_cc $energy  >> curve.amber.c2h4
#
# loop over n
#
done
 
#
# Cleaning...
#
mkdir Logs
mv c2h4.*.inp c2h4.*.out Logs
rm RUN* *restart*
 
  • Use a fitting gnuplot script to verify the value of the parameters. What has changed with respect to the C2H2 case? Why?
Assignment:
  1. Report the values of the bond length and force constants in all cases.
  2. Discuss these values in view of the kind of molecules we are simulating.
exercises/2015_ethz_mmm/c2h2_bond_energy.1425421546.txt.gz · Last modified: 2015/03/03 22:25 by dpasserone