no way to compare when less than two revisions

Differences

This shows you the differences between two versions of the page.

@@ Line 1: / Line 1: @@
+====== 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 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 (topology file)
+      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 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 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 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>