exercises:2017_ethz_mmm:c2h2_bond_energy
Differences
This shows you the differences between two versions of the page.
— | exercises:2017_ethz_mmm:c2h2_bond_energy [2020/08/21 10:15] (current) – created - external edit 127.0.0.1 | ||
---|---|---|---|
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 important> | ||
+ | |||
+ | < | ||
+ | you@eulerX ~$ module load 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:// | ||
+ | you@eulerX ~$ unzip exercises: | ||
+ | </ | ||
+ | |||
+ | <note tip> | ||
+ | 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 " | ||
+ | < | ||
+ | you@eulerX ~$ cd exercise_2.1/ | ||
+ | </ | ||
+ | The input file structure of the template is the following: | ||
+ | |||
+ | <code - c2h2.inp.templ> | ||
+ | & | ||
+ | METHOD FIST ! Using Molecular Mechanics | ||
+ | &MM | ||
+ | & | ||
+ | parm_file_name c2h2.top | ||
+ | ! This file contains force field parameters (topology file) | ||
+ | parmtype AMBER | ||
+ | & | ||
+ | EMAX_SPLINE 15.0 ! Specify the maximum value of the potential up to which splines will be constructed | ||
+ | &END | ||
+ | &END FORCEFIELD | ||
+ | & | ||
+ | & | ||
+ | EWALD_TYPE SPME ! Smooth particle mesh using beta-Euler splines | ||
+ | ALPHA .36 | ||
+ | GMAX 128 ! Number of grid points | ||
+ | &END EWALD | ||
+ | &END POISSON | ||
+ | & | ||
+ | & | ||
+ | &END | ||
+ | & | ||
+ | &END | ||
+ | &END | ||
+ | &END MM | ||
+ | & | ||
+ | & | ||
+ | ABC 50 50 50 ! Lengths of the cell vectors A, B, and C | ||
+ | &END CELL | ||
+ | & | ||
+ | H 0 0 0 | ||
+ | C | ||
+ | C | ||
+ | H | ||
+ | & | ||
+ | & | ||
+ | CHARGE_BETA | ||
+ | CONNECTIVITY AMBER ! Use AMBER topology file for reading connectivity | ||
+ | CONN_FILE_NAME c2h2.top | ||
+ | ! This file contains connectivity data | ||
+ | &END TOPOLOGY | ||
+ | |||
+ | & | ||
+ | AMBER_INFO | ||
+ | &END | ||
+ | &END | ||
+ | &END SUBSYS | ||
+ | &END FORCE_EVAL | ||
+ | |||
+ | & | ||
+ | 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 | ||
+ | &END GLOBAL | ||
+ | </ | ||
+ | |||
+ | |||
+ | 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 | ||
+ | . / | ||
+ | |||
+ | |||
+ | # | ||
+ | # 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 ' | ||
+ | echo $dist_cc $energy | ||
+ | # | ||
+ | # 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 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 | ||
+ | </ | ||
+ | |||
+ | <code - fit.gnu> | ||
+ | # | ||
+ | f(x)=k*(x-x0)*(x-x0)+e0 | ||
+ | delta = 0.002 # | ||
+ | e0=0.007720330647107 | ||
+ | k = 500 # initial value of k | ||
+ | x0=1 # initial value of x0 | ||
+ | #a0 = 1.3327436840 | ||
+ | set xlabel 'd (Angstrom)' | ||
+ | set ylabel 'E (kcal/ | ||
+ | set title "c2h2 triple bond" | ||
+ | |||
+ | fit f(x) ' | ||
+ | # 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 ' | ||
+ | # ^ 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" | ||
+ | |||
+ | 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: | ||
+ | |||
+ | <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 | ||
+ | </ | ||
+ | |||
+ | * 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 | ||
+ | #. ~/ | ||
+ | . / | ||
+ | |||
+ | |||
+ | # | ||
+ | # 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 ' | ||
+ | echo $dist_cc $energy | ||
+ | # | ||
+ | # 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> | ||
+ | - 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. | ||
+ | </ |