Differences

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

--- exercise:uv [2014/05/16 02:56] – created dpasserone
+++ exercises:2014_ethz_mmm:uv [2020/08/21 10:15] (current) – external edit 127.0.0.1
@@ Line 12: / Line 12: @@
 </code>
-Then source your profile file:
+Then source your profile file, as well as loading the modules, and copying this configuration file:
 <code>
 . ~/.bash_profile
+module load intel/12.1.2 open_mpi/1.6.5 python vmd
+cp ~danielep/.nwchemrc $HOME
 </code>
-Now you are able to run the **nwchem** code. Theo
+Now you are able to run the **nwchem** code.
+Download all files from the media manager: {{exercise_11.1.tar.gz|}}.
+===== Calculation of the spectrum using linear response TDDFT =====
-We will show how a simple change in the termination (1 vs. 2 Hydrogens) changes the state structure completely.
+We first run a calculation using linear response TDDFT. The code computes the ground state, then computes all excitations with their dipole oscillator strength, keeping into account open-shell configurations. We use the hybrid functional PBE0 here. The input file has a simple structure and includes coordinates, basis set specification and level of DFT.
-{{ :exercise:2rib.jpg?direct&300 |}}
+<code>
-<note tip>
+# This tests CIS, TDHF, TDDFT functionality at once
-You should run these calculations on 16 nodes with ''bsub -n 16''.
+# by using a hybrid LDA, GGA, HF functional for
-Copy, as usual, the files from the directory **/cluster/home03/matl/danielep/LECTURE10/EXERCISE_10.2** (and later here on the media manager).
+# spin restricted reference with symmetry on.
-</note>
-===== 1. Task: Running the job and looking at the orbitals =====
+start tddft_h2o
-This time we will not optimize the structure. With an ENERGY run, we run with ** cp2k ** the job 1h.1.5.inp and 2h.1.5.inp, meaning that there are here 1.5 units of the original molecule in the gas phase.
-The interesting part of the code is in the ** &PRINT ** section of ** &DFT **:
+echo
+title "TDDFT H2O PBE0/6-31G**"
+geometry
+O     0.00000000     0.00000000     0.12982363
+H     0.75933475     0.00000000    -0.46621158
+H    -0.75933475     0.00000000    -0.46621158
+end
+basis
+O library 6-31G**
+H library 6-31G**
+end
+dft
+xc pbe0
+odft
+end
+dplot
+   TITLE h2o
+   LimitXYZ
+   -2.0 2.0 25
+   -2.0 2.0 25
+   -2.0 2.0 25
+   spin total
+   orbitals view; 1; 7
+   gaussian
+end
+tddft
+ nroots 20
+end
+task tddft energy
+task dplot
+</code>
+The command is
 <code>
-    &PRINT
+bsub  -n 2 -o tddft_h2o_uhf.out  mpirun nwchem tddft_h2o_uhf.nw
-    &STM
+</code>
-      BIAS -2.0 -1.0 1.0 2.0
-      TH_TORB S
-    &END STM
-    &MO_CUBES
-       NHOMO 10
-       NLUMO 10
-       STRIDE 2 2 2
-       WRITE_CUBE T
-    &END
-    &V_HARTREE_CUBE
-    &END
-      &MO
+Once the job has started, you can monitor the output file by the command (the file is still not present in your directory):
-       FILENAME EIG
-       ADD_LAST NUMERIC
+<code>
-       EIGENVALUES
+bpeek -f [jobid]
-       OCCUPATION_NUMBERS
-      &END
-    &END
 </code>
-There will be an output file with the energy levels and their occupation. The last one can be easily found...
+where **jobid** is the job number (see bjobs) and is necessary only if you are running more than one job.
+In the output file (**  tddft_h2o_uhf.out ** ) you can find all orbital energies.
+You can also plot the orbitals with vmd. There are lumo.cube and homo.cube files generated previously with the input **uhf.nw**.
+To visualize them,
-<note tip>
+<code>
-Hitting ** ls -ltr ** will allow you to see on the last lines of the screen the most recent files.
+vmd -e homo.vmd
-</note>
+vmd -e lumo.vmd
+</code>
+Then, the excitation spectrum can be visualized using
+<code>
+python ./nwchem_tddft_spectrum.py -o tddft_h2o_uhf.out
+</code>
 <note important>
-  - Draw the energy level diagram for the two molecules. What is the energy gap in the two cases? What are the differences?
+  - List in a table the orbital energies for this system. Note that alpha and beta orbitals are listed, but they are degenerate in this case (alpha=beta). Search for the string "Occ." just after the title "Final Molecular Orbital Analysis". Note the character of the orbital (px, py, pz...) and the sign of the LUMO coefficients (comment).
-  - Look with vmd at the cube files corresponding to the most interesting levels (close to Fermi...). Comment on the distribution of the states.
+  - Visualize homo and lumo with vmd.
+  - The excitation spectrum corresponds to transitions between occupied and unoccupied states. Look for this information in the file, and compare it with the peaks in the plot.
 </note>
-===== 2. Task: Producing a simple STM image =====
+===== Resonant ultraviolet excitation of water =====
-The section ** &STM ** shown above produces STM images at different bias (feel free to change), meaning, using the Tersoff-Hamann approximation, it integrates all the density of states with energies between Fermi energy and the Bias potential: this energy interval is involved in the tunnel current.
+In this second part, we compute the time-dependent electron response to a quasi-monochromatic laser pulse tuned to a particular transition.
-The ** *STM*cube files are 3D maps of the integrated density of states. Imagine that we have a microscope with a feedback that can keep constant current between tip and sample, by changing the height of the tip on the surface. Since the current is proportional to the density of states, we move the tip on ** isosurfaces ** of our cubefile.
+The spectrum obtained with linear response TDDFT can be also calculated by exciting the system through a laser pulse with a specific polarization along x, y, or z.
-The program ** stm ** (in the same working dir) allows to extract a 2D map of the height of a given isosurface.
+We will use the results of a calculation described [[http://www.nwchem-sw.org/index.php/Release62:RT-TDDFT|here]]. The total spectrum  (6-31G * */PBE0 gas-phase water) is the same as the one we calculated before. First, we consider the absorption spectrum (computed previously) but plotted for the three polarizations (x,y,z) rather then as a sum. The details are given in the previously cited link.
+{{ 730px-rt_tddft_h2o_resonant_spec_field.png?direct&300 |}}
+Say we are interested in the excitation near 10 eV. We can clearly see this is a z-polarized transition (green on curve). To selectively excite this we could use a continuous wave E-field, which has a delta-function, i.e., single frequency, bandwidth but since we are doing finite simulations we need a suitable envelope. The broader the envelope in time the narrower the excitation in frequency domain, but of course long simulations become costly so we need to put some thought into the choice of our envelope. In this case the peak of interest is spectrally isolated from other z-polarized peaks, so this is quite straightforward. The procedure is outlined below, and the corresponding frequency extent of the pulse is shown on the absorption figure in orange. Note that it only covers one excitation, i.e., the field selectively excites one mode. The full input deck is ** h2o_resonant.nw **.
+The relevant code section is:
 <code>
-Run the program in the following way:
+rt_tddft
-$ module load boost/1.54.0
+  tmax 1000.0
-$ module load mkl
+  dt 0.2
-$ stm -c 2h*STM*.cube --isovalues 1E-5 > stm.out
-</code>
-The resulting .igor files contain the z profile (in bohr) and may for example be plotted by gnuplot:
+  field "driver"
+    type gaussian
+    polarization z
+    frequency 0.3768  # = 10.25 eV
+    center 393.3
+    width 64.8
+    max 0.0001
+  end
+  excite "system" with "driver"
+ end
+task dft rt_tddft
+</code>
+Run now ** h2_resonant.nw ** on 4 cores:
 <code>
-gnuplot
+bsub -n 4 -o resonant.out  mpirun nwchem resonant.nw
-set pm3d map
+</code>
-set size square
-set xrange [......
-set yrange [.....
-splot "mystm.igor" matrix using 2:1:3
+The run (follow with bpeek) will apply a field for a limited amount of time. This field will excite the system into a superposition of the ground state and the one excited state, which manifests as monochromatic oscillations. After the field has passed the dipole oscillations continue forever as there is no damping in the system.
+There are 5000 MD steps. It will take about 10 minutes. At the end, cubefiles for the density at each timestep will be generated.
+You can visualize the animation using vmd, using a script that cleans a bit (moving the cubes into another directory)
+<code>
+./createlist
+vmd -e animate.cube.vmd
+</code>
-</code>
+What you will see is the electron density difference between the initial state and an instant along the trajectory.
-Where instead of "mystm" you use an appropriate filename.
 <note important>
-  - In the output file of cp2k, the program tells you how many states have contributed to each STM image. Discuss the images that you see in the two cases.
+  - Plot from the output file the applied field: **grep -i Applied resonant.out | grep alpha > appl**
-  - What makes the 1h* case particular with respect to the 2h*?
+  - Plot the z component of the induced dipole moment: **grep ipole resonant.out > dipole**
-  - Change the isosurface and look at the z-range. Discuss the changes in the range.
+  - Explain what you see in the vmd representation based on what you see on the previous plot
-  - Would you define the differences between 1h and 2h in the STM images as more of structural origin or electronic origin?
 </note>