User Tools

Site Tools


exercises:2017_ethz_mmm:c2h2_bond_energy

Differences

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

Link to this comparison view

exercises:2017_ethz_mmm:c2h2_bond_energy [2017/02/22 10:01]
exercises:2017_ethz_mmm:c2h2_bond_energy [2017/02/22 10:01] (current)
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>​
exercises/2017_ethz_mmm/c2h2_bond_energy.txt ยท Last modified: 2017/02/22 10:01 (external edit)