====== 38 atom Lennard-Jones cluster ====== {{:exercises:2017_ethz_mmm:lj38bs.jpg?400|}} (picture by Luke Abraham) All files of this exercise be downloaded directly from the wiki: {{exercise_1.1.zip|}} Download the 1.1 exercise into your **EXERCISES** folder and unzip it. max@qmobile:~$ cd ; cd EXERCISES max@qmobile:~$ wget http://www.cp2k.org/_media/exercises:2018_ethz_mmm:exercise_1.1.zip max@qmobile:~$ unzip exercises:2018_ethz_mmm:exercise_1.1.zip max@qmobile:~$ cd exercise_1.1 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: [[doi>10.1063/1.478595]] The command to run cp2k is the following (with a generic **file.inp** input file): max@qmobile:~$ cp2k.ssmp -i file.inp -o file.out ===== 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 input file structure of the template is the following: &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 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 K_b to make units work...). 1 Kelvin*K_b=3.2E-6 Hartree. Using this conversion factor you can transform the epsilon value into Hartree, and the total energy can be expressed in units of epsilon. **The sigma value is in Angstrom.** - randomize the coordinate files **fcc.xyz** (which represents the "cubic" structure) m_xyzrand 1.0 < fcc.xyz > fcc_rand.xyzDo the same with **ico.xyz** which represents the icosahedral structure. You can look at all files with **vmd**. - extract the q4 order parameter from **fcc.xyz** and from **fcc_rand.xyz** and compare the values. - python stein.py file.xyz You will be asked the cutoff radius for the neighbors, it is **1.391** in sigma units. **You should input it in Angstrom**. You will also be asked **"value of l"** This means the symmetry of the order parameter, which is **l=4** in this case. - before running the simulation, copy the input coordinate file into in.xyz cp fcc_rand.xyz in.xyz - Before running cp2k, check if the file **OPT-pos-1.xyz** is already present from a previous run. In that case remove or delete it accordingly. It contains the trajectory of the optimization. - run cp2k cp2k.ssmp -i geo_opt.inp | tee geo_opt.out (to see the output on the screen as well), or **AS AN ALTERNATIVE** cp2k.ssmp -i geo_opt.inp > geo_opt.out (to retain the output in the geo_opt.out file only) - in the output file, grep the final energy grep "ENERGY|“ geo_opt.out and transform it in the unit of the paper (epsilon units) - Open vmd and play with the optimization trajectory vmd OPT-pos-1.xyz (ask the teacher) - apply the script **myq4** to the optimization trajectory: this generates a list of q4 and energies for the whole trajectory. ./myq4 OPT-pos-1.xyz > fcc.ene.q4 - plot q4 and energies with **gnuplot** (ask the teacher) - have a look at the myq4 script nano myq4 - repeat for the ico.xyz starting point, don't forget to first copy/remove the files appropriately. For example: mkdir FCC ; mv OPT* FCC ; mv geo_opt.out FCC - Run the bash script ./curveLook inside, and try to understand what you get. - create a FCC_OUT subdirectory (**mkdir FCC_OUT ; cd FCC_OUT**) and copy there the files you want to keep; then go back one dir (**cd ..**), delete all the OPT* files (**rm OPT* **) and repeat the exercise with ico.xyz Assignment: - Report the energy of the minima, compare it with the ones of the initial configurations. - After converting the energy into "epsilon" units, estimate the number of bonds in the cluster, assuming a pairwise interaction. - 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.