This is an old revision of the document!
C2H2 and C2H4 bond energy
you@eulerX ~$ ./update_functionlibrary
If you don't have the file “update_functionlibrary” please copy it using the following command:
you@eulerX ~$ cp ~yakutova/update_functionlibrary ~/
or download it from here update_functionlibrary and put it into ~/Scripts/ folder
you@eulerX ~$ module load new cp2k
and to submit the job chain:
you@eulerX ~$ bsub < c2h2.a.chain
Create a new directory for this exercise:
you@eulerX ~$ mkdir mmm_exercise_2.1 you@eulerX ~$ cd mmm_exercise_2.1
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.
Please copy them to your directory. The input file structure of the template is the following:
- c2h2.amber.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.amber.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@brutusX mmm_exercise_2.1$ module load new cp2k you@brutusX mmm_exercise_2.1$ bsub < c2h2.a.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@brutusX mmm_exercise_2.1$ gnuplot fit.triple.gnu
- fit.triple.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.
- The input c2h4.amber.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.a.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?
<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