Differences

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

--- exercises:2017_ethz_mmm:lennard_jones_cluster [2017/02/22 06:14] – external edit 127.0.0.1
+++ exercises:2017_ethz_mmm:lennard_jones_cluster [2020/08/21 10:15] (current) – external edit 127.0.0.1
@@ Line 1: / Line 1: @@
 ====== 38 atom Lennard-Jones cluster ======
+{{:exercises:2017_ethz_mmm:lj38bs.jpg?400|}} (picture by Luke Abraham)
 <note warning>
 TO USE THE FUNCTION LIBRARY (VERSION UP TO DATE) IN THE INTERACTIVE SHELL:
@@ Line 16: / Line 18: @@
 </code>
-**and to submit the job:**
+**and to submit the job (note: since all the examples of this week are ultrafast, we will run them interactively, and NOT on a compute node. This is not the normal procedure for the next lectures).**
 <code>
-you@eulerX ~$ bsub < jobname
+you@eulerX ~$ cp2k.popt -i file.inp -o file.out
 </code>
 </note>
@@ Line 31: / Line 33: @@
 <note tip>
-All files of this exercise (**input and scripts are all commented**) can be downloaded from the wiki: {{exercise_1.1.zip|}}
+All files of this exercise be downloaded from the wiki: {{exercise_1.1.zip|}}
 </note>
-In this exercise you will test the Lennard-Jones potential. Login to euler using your nethz credentials.
+In this exercise you will test the Lennard-Jones potential. In particular, we will focus on the system described in the following paper about the energy landscape of the 38 atom Lennard-Jones cluster:
-Then go to the directory "exercise_1.1/38/".
+<note tip>[[doi>10.1063/1.478595]]
+</note>
+Login to euler using your nethz credentials.
+Then go to the directory "exercise_1.1".
 <code>
-you@eulerX ~$ cd exercise_1.1/38/
+you@eulerX ~$ cd exercise_1.1
-</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>
+===== Geometry optimization  =====
+In this first part you will perform a simple energy optimization, to find the two lowest lying minima in the potential energy surface.
-The script to modify the C-C distance in the input and generate the curve is the following:
+The input file structure of the template 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 - geo_opt.inp>
+&GLOBAL
+ FLUSH_SHOULD_FLUSH
+ PRINT_LEVEL low
+ PROJECT geo_opt_bfgs
+ RUN_TYPE geo_opt
+ WALLTIME 600
+&END GLOBAL
+&MOTION
+ &GEO_OPT
+  OPTIMIZER BFGS
+  MAX_ITER  200
+  MAX_DR    0.001
+  RMS_DR    0.0003
+  MAX_FORCE 0.0001
+  RMS_FORCE 0.00003
+  &BFGS
+   USE_MODEL_HESSIAN yes
+  &END BFGS
+ &END GEO_OPT
+ &PRINT
+  &TRAJECTORY on
+   FORMAT xyz
+   &EACH
+    GEO_OPT 1
+   &END EACH
+  &END TRAJECTORY
+ &END PRINT
+&END MOTION
+&FORCE_EVAL
+ METHOD Fist
+ STRESS_TENSOR ANALYTICAL
+ &MM
+    &FORCEFIELD
+      &CHARGE
+        ATOM Ar
+        CHARGE 0.0
+      &END
+      &NONBONDED
+        &LENNARD-JONES
+          atoms Ar Ar
+          EPSILON 119.8
+          SIGMA 3.405
+          RCUT 8.4
+        &END LENNARD-JONES
+      &END NONBONDED
+      &CHARGE
+        ATOM Kr
+        CHARGE 0.0
+      &END CHARGE
+    &END FORCEFIELD
+  &POISSON
+   PERIODIC NONE
+   &EWALD
+    EWALD_TYPE none
+   &END EWALD
+  &END POISSON
+  &PRINT
+   &FF_INFO OFF
+    SPLINE_DATA
+    SPLINE_INFO
+   &END FF_INFO
+  &END PRINT
+ &END MM
+ &PRINT
+  &FORCES off
+  &END FORCES
+  &GRID_INFORMATION
+  &END GRID_INFORMATION
+  &PROGRAM_RUN_INFO
+   &EACH
+    GEO_OPT 1
+   &END EACH
+  &END PROGRAM_RUN_INFO
+  &STRESS_TENSOR
+   &EACH
+    GEO_OPT 1
+   &END EACH
+  &END STRESS_TENSOR
+ &END PRINT
+ &SUBSYS
+  &CELL
+   A      100 0 0
+   B      0   100 0
+   C      0 0 100
+   PERIODIC NONE
+  &END CELL
+  &TOPOLOGY
+      COORD_FILE_NAME in.xyz
+      COORDINATE xyz
+  &END
+  &PRINT
+   &CELL
+   &END CELL
+   &KINDS
+   &END KINDS
+   &MOLECULES OFF
+   &END MOLECULES
+   &SYMMETRY
+   &END SYMMETRY
+  &END PRINT
+ &END SUBSYS
+&END FORCE_EVAL
 </code>
+<note important>NOTE ON THE UNITS: CP2K USES SO CALLED "atomic units". Meaning that the resulting energies are expressed in Hartree,
+**1 Hartree=27.2114 eV**.
+In the input file, the epsilon value (depth of the well) is expressed in KT units, namely, in "temperature" units (there is a Boltzmann constant to make units work...). **The sigma value is in Angstrom.**
+</note>
+<note tip>
+  - load the module with the special m_* bash functions and initialize the module: <code>module load courses mmm ; mmm-init </code>
+  - randomize the coordinate files **fcc.xyz**  <code>m_xyzrand 1.0 < fcc.xyz > fcc_rand.xyz</code> . Do the same with ico.xyz
+  - extract the q4 order parameter from **fcc.xyz** and from **fcc_rand.xyz** and compare the values.<code>module load new gcc/4.8.2 python/2.7.12
+python stein.py file.xyz </code>. You will be asked the cutoff radius for the neighbors, it is **1.391** in sigma units. **You should input it in Angstrom**.
+  - before running the simulation, copy the input coordinate file into in.xyz <code>cp fcc_rand.xyz in.xyz</code>
+  - run cp2k <code>module load cp2k</code>(this has only to be done once)<code>cp2k.popt -i geo_opt.inp -o geo_opt.out </code>
+  - in the output file, note the final energy, **transform it in the unit of the paper (epsilon units)**
+  - load vmd module and play with the optimization trajectory <code>vmd OPT-pos-1.xyz</code> (ask the teacher)
+  - apply the script **myq4** to the optimization trajectory: this generates a list of q4 and energies for the whole trajectory. <code>./myq4 OPT-pos-1.xyz > fcc.ene.q4</code>
+  - plot q4 and energies with **gnuplot** (ask the teacher)
+  - have a look at the myq4 script <code>nano myq4</code>
+  - repeat for the ico.xyz starting point, don't forget to first copy/remove the files appropriately. For example: <code>mkdir FCC ; mv OPT* FCC ; mv geo_opt.out FCC</code>
+  - finally, run the bash script <code>./curve</code>. Look inside, and try to understand what you get.
+</note>
-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.
+  - Report the energy of the minima, compare it with the ones of the initial configurations.
-  - Discuss these values in view of the kind of molecules we are simulating.
+  - Plot q4 vs. energy and q4 vs. optimization steps, for the two cases. Discuss the results. Are the minima in two separate basins?
+  - Report the value of the order parameter of the minumum, and discuss what you see
+  - Use "gnuplot" to make the output of "./curve" understandable, discuss the results.
 </note>