NWChem computer exercise

Quantumchemistry, NWI-MOL406, Computer exercises, week 5 (13-Dec-2012)

Calculate and animate vibrations of Furan

Follow the installation instructions from last week
Before you start ECCE define your favorite editor. The default is /usr/bin/vim. If you want to use /usr/bin/gedit instead (in tcsh):
```
setenv EDITOR /usr/bin/gedit
```
Note: to find out where a program is use:
```
which gedit
```
Now start ecce:
```
ecce &
```
We will use the "organize" to prepare the input for NWChem, to start the calculation, and to view the results:
After logging into the ECCE Gateway, click on "organizer"
Click in the menu bar on "HELP" and choose "Help on organizer", and read it
The first goal is to animate vibrations in Furan, calculated at the HF/STO-3G level

Here are some hints to navigate through the menus:

organizer
  right click on user
    new NWchem calculation
      editor
        builder
          structure-library
            SimpleStructure
               Simple_heterocyclic_compounds
                 Furan
Work area:
  click to Add Structure
    save
      File: quit

ECCE NWChem Editor
  Basis Set
    Quick Basis Menu
      STO-3G (29 functions, 87 primitives)

NWCHem Settings
   Theory RHF (default)
   Runtype: Vibration

Details: have a quick look here

Final Edit: click to see whether /usr/bin/gedit of /usr/bin/vi will start

Launch
 Before you can start you have to tell ECCE which computer to use:
 Click on "Job" in the menu bar and select:
   Register Machines
     Machine: localhost   (you can also use hg206pc12.science.ru.nl, or whatever PC you are using)
       check "ssh"
         Applications:
         NWCHem: /vol/xpcursus/QuantumChemistry/bin/nwchem
       Add/Change
       Close

Go back to ECCE Launcher:
  Machine: localhost
     Username:  xxx  Password: (the password to login into Ubuntu)
     Remote Shell: ssh
     click on "Open" to make sure "ssh" can login with the password
  Run directory: /scratch/gerritg (or whatever your username is)

Launch
  last message should be: Successfully submitted job

Close ECCE NWCHem Editor

Already while the job runs you can use this button to follow the
calculation:
ECCE Organizer: Run Mgmt
      - shell in calculation run directory (should start a terminal)
      - View Run log
      - View Output file

ECCE Organizer: Click Viewer
  Vibrational frequencies
    Select Animation

You should now be able to select a vibrational frequency and animate the corresponding vibrational mode. It is also possible to look at the table of frequencies.
Right click in the spectrum and select "Clone current plot" to open the spectrum in a larger window.

Calculate and animate vibrations of Furan: correcting the mistake

Before calculating vibrational frequencies we should have optimized the geometry. We will do this in a new NWChem run, so we can see how much difference it makes:

In the ECCE Organizer project tree right click to copy, paste, and rename your last run.
Start the calculation editor, and click on the blue reset button in the lower left corner.
Use the "Final edit" button to change the input: instead of the last two lines:
```
task scf energy
task scf freq
```
use
```
task scf optimize
task scf freq
```
Now run the calculation again. Since the geometry optimization takes a little longer, you may see the "tail -f on output" option in the "Run Mgmt" menu. With this option you can look at the NWChem output while it is generated.
When it is is done, use the "ECCE Viewer" to look at the new frequencies. Use two Viewer windows to put the old an new frequencies next to each other. Does the geometry optimization make the frequencies higher or lower?

Anthracene in D2h symmetry

NWChem/ECCE does not always recognize the proper symmetry. The builder has a "symmetry tool" to set the proper symmetry by hand.

The D2h character table
Start a new calculation of anthracene

In the builder/structure library it can be found in

builder
  structure-library
    SimpleStructure
      Simple_arenes
        Anthracene
Add structure in work area

Now switch on "Tools/symmetry"
Select D2h and click find
Also click "Clean symmetry", otherwise the symmetry may be reduced to, e.g., C2v in the calculation
Do a SCF calculation in the STO-3G basis
Plot the HOMO

Recalculate this orbital with the Huckel method and compare the result

Use the "Options" to switch off the Hydrogen atoms
Switch on "atom labels", "elements with number"
Start Matlab and do a Huckel calculation for anthracene:

First make a table of neighbors:

tab = [ ...
 1   2
 2   3
 3   4
etc...
]

Determine the number of bonds (should be 16):
```
m = size(tab,1)
```

Now set up the Hamiltonian matrix:

n   = 14  % dimension of Hamiltonian
H   = zeros(n, n);
for i=1:m
  H(tab(i,1), tab(i,2)) = 1;
end
H = H + H'    % make it Hermitian

There is a quicker way to setup the Hamiltonian. After you enter the table use "sparse":
```
H = sparse(tab(:,1), tab(:,2), 1)
H = H+H'
H = full(H)
```
Now calculate eigenvalues and eigenvectors
```
[V,E] = eig(H);
```
The eigenvalues should be on the diagonal of E, so it is convenient to put them in a vector
```
e = diag(E)
```
Now think long and hard which column of V is the HOMO
```
k = 8
```
Make convenient table:
```
[ (1:n)' V(:,k)]
```
Use "View coefficient" and compare the result with the HF/STO-3G calculation
Copy the coefficients of the Px orbitals in a column vector, and determine the overlap with all the Huckel MOs:
```
mo = [0.29
      0.23
      etc]

S = V'*mo
```
Is is also possible to "export" the coefficients as a "comma separated values" (.csv) file (right-click on table). To import the file in Matlab use "File", "Open" and select "Exclude columns with", "Non-numeric cells".
Import the file into the matrix A (80x2)
To select the non-zero entries use:
```
ind = find(A(:,2)~=0)
mo   = A(ind, 2)
```
Or even shorter:
```
ind = find(A(:,2))
mo  = A(ind,2)
```
Use the Viewer to find an orbital of B3g symmetry
In Matlab define all SALCs of B3g in a matrix, where each column corresponds to a SALC:
```
D = [ ...
1 0 0
0 1 0
0 0 1
etc
]
```

Normalize each column:

l = size(D,2);
for i=1:l
  D(:,i) = D(:,i)/norm(D(:,i));
end

Check the orthonormality of the basis:
```
S = D'*D
```
Set up the Hamiltonian in the symmetry basis:
```
Hb3g = D'*H*D
```

Solve Huckel in symmetry basis:

[Vb3g, Eb3g] = eig(Hb3g);
eb3g = diag(Eb3d)

Transform the eigen vector back to the origional basis:
```
Wb3g = D * Vb3g
```
Compare with HF/STO3g results.

Do not forget de-install before you log-off.

Back to [ Quantum Chemistry ]

Last updated: 13-Dec-2012, by Gerrit C. Groenenboom.