CP2K Open Source Molecular Dynamics

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== C2H2 and C2H4 bond energy ======
<note warning>
TO USE THE FUNCTION LIBRARY (VERSION UP TO DATE) IN THE INTERACTIVE SHELL:

you@eulerX ~$ module load courses mmm vmd
 
you@eulerX ~$ mmm-init
</note>



<note important> **REMEMBER: this is the command to load the  module for the cp2k program:**

<code>
you@eulerX ~$ module load new cp2k
</code>

**and to submit the job:**

<code>
you@eulerX ~$ bsub < jobname
</code> 
</note> 

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

<code>
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
</code>

<note tip>
All files of this exercise (**input and scripts are all commented**) can be downloaded from the wiki: {{exercise_2.1.zip|}} 
</note>

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/". 
<code>
you@eulerX ~$ cd exercise_2.1/c2h2/
</code>
The input file structure of the template is the following:

<code - 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                                                                                                                                                                                
</code>


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

<code bash 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: m_change m_replace m_getcolumn m_sum m_list
#
module load new cp2k
. /cluster/apps/courses/mmm/m_functions.bash


#
# deletes old output file
#
rm curve.amber.c2h2

#
# loop over all values of n
#

for n in `m_list -5 5`
do
  dist_cc=`m_change $dist_cc0 $n $delta`
  c2=`m_sum $dist_hc $dist_cc`
  h2=`m_sum $c2 $dist_hc `
  #
  # replaces coordinates in the input
  #
  m_replace _C1_ $dist_hc < c2h2.inp.templ | m_replace _C2_ $c2  | m_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=` m_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*



</code>


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

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



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

<code>
you@eulerX c2h2$ gnuplot fit.gnu
</code>

<code - 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
</code>

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.
<code>
you@eulerX c2h2$ cd ../c2h4
</code>

    * The new input c2h4.inp.templ is similar to the previous one, only with different project name (C2H4) and coordinates:

<code coo.ch4>
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 
</code>

   * The script to modify the bond lengths changes slightly as well

<code bash 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: m_change m_replace m_getcolumn m_sum m_list
#
module load new cp2k
#. ~/Scripts/myfunctions.bash
. /cluster/apps/courses/mmm/m_functions.bash


#
# deletes old output file
#
rm curve.amber.c2h4

#
# loop over all values of n
#

for n in `m_list -5 5`
do
  dist_cc=`m_change $dist_cc0 $n $delta`
  c2=`m_sum $dist_hcx $dist_cc`
  h3=`m_sum $c2 $dist_hcx `
  #
  # replaces coordinates in the input
  #
  m_replace _C1_ $dist_hcx < c2h4.inp.templ | m_replace _C2_ $c2  | m_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=` m_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*

                                                                                                                                                  
</code>

  * Use a fitting gnuplot script to verify the value of the parameters. What has changed with respect to the C2H2 case? Why?


<note tip>Assignment: 
  - Report the values of the bond length and force constants in all cases. 
  - Discuss these values in view of the kind of molecules we are simulating.
</note>