646. NWChem 6.6 on debian stretch (9)

I've got plenty of posts on how to compile nwchem on this blog. However, as I'm still running debian jessie on my work computers and cluster I have failed to appreciate quite how different things work in debian stretch.

This includes the compilation of ECCE (although owing to updates meant to rectify that at https://github.com/FriendsofECCE/ECCE this isn't necessarily something you'd notice) and NWChem.

Anyway, here's how to build nwchem with openmpi support under debian stretch:

---

1. Download source code and patches

sudo mkdir /opt/nwchem -p
sudo chown $USER:$USER /opt/nwchem
cd /opt/nwchem
mkdir patches6.6
cd patches6.6
for i in Tddft_mxvec20.patch.gz Tools_lib64.patch.gz Config_libs66.patch.gz Cosmo_meminit.patch.gz Sym_abelian.patch.gz Xccvs98.patch.gz Dplot_tolrho.patch.gz Driver_smalleig.patch.gz Ga_argv.patch.gz Raman_displ.patch.gz Ga_defs.patch.gz Zgesvd.patch.gz Cosmo_dftprint.patch.gz Txs_gcc6.patch.gz Gcc6_optfix.patch.gz Util_gnumakefile.patch.gz Util_getppn.patch.gz Gcc6_macs_optfix.patch.gz Xatom_vdw.patch.gz Notdir_fc.patch.gz Hfmke.patch.gz Cdfti2nw66.patch.gz ; do wget http://www.nwchem-sw.org/download.php?f=$i ; done
gunzip *.gz
cat *.patch > consolidated.patch
cd ..
wget http://www.nwchem-sw.org/download.php?f=Nwchem-6.6.revision27746-src.2015-10-20.tar.gz -O Nwchem-6.6.revision27746-src.2015-10-20.tar.gz
tar xvf Nwchem-6.6.revision27746-src.2015-10-20.tar.gz

Apply patches:

cd /opt/nwchem/nwchem-6.6/
cp /opt/nwchem/patches6.6/consolidated.patch .
patch -p0 < consolidated.patch

In addition, you can do the following to print more information so that GabEdit properly shows MOs:

cd nwchem-6.6/
grep -rl "do klo = 0, min(n-1,9), 2"|xargs -I {} sed -i 's/do klo = 0, min(n-1,9), 2/do klo = 0, min(n-1,199), 2/g' {}
grep -rl " 0.15d0,"|grep -v "atomdata"|xargs -I {} sed -i 's/0.15d0/0.01d0/g' {}

---

2. Install dependencies and compile.
 How exactly it is done depends on your compiler and math libs. I've used red for examples with INTEL mkl, blue for AMD acml and orange for the free openblas. ifort is shown in purple and gfortran like this.You should only use one compiler + one math lib.
Hint: also make sure your ld.conf.so.d files contain paths to your math libs.

cd /opt/nwchem/nwchem-6.6
sudo apt-get install build-essential gfortran python2.7-dev libopenmpi-dev openmpi-bin

export NWCHEM_TOP=`pwd`
export LARGE_FILES=TRUE
export TCGRSH=/usr/bin/ssh
export NWCHEM_TOP=`pwd`
export NWCHEM_TARGET=LINUX64
export NWCHEM_MODULES="all"
export PYTHONVERSION=2.7
export PYTHONHOME=/usr

export BLASOPT="-L/opt/intel/composer_xe_2013.4.183/mkl/lib/intel64/ -lmkl_core -lmkl_sequential -lmkl_intel_ilp64"
export LIBRARY_PATH="$LIBRARY_PATH:/opt/intel/composer_xe_2013.4.183/mkl/lib/intel64/"

export BLASOPT="-L/opt/acml/acml5.3.1/gfortran64_int64/lib -lacml"
export LIBRARY_PATH="$LIBRARY_PATH:/opt/acml/acml5.3.1/gfortran64_int64/lib"

export BLASOPT="-L/opt/openblas/lib -lopenblas"
export LIBRARY_PATH="$LIBRARY_PATH:/opt/openblas/lib"


export USE_MPI=y
export USE_MPIF=y
export USE_MPIF4=y
export ARMCI_NETWORK=SOCKETS

cd $NWCHEM_TOP/src

make clean
make nwchem_config
make FC=ifort 1> make.log 2>make.err
make FC=gfortran 1> make.log 2>make.err

cd $NWCHEM_TOP/contrib
export FC=ifort
export FC=gfortran
./getmem.nwchem

---

3. Run
mpirun -n 8 /opt/nwchem-6.6/bin/LINUX64/nwchem input.nw > output.nw

630. Making MO figures in gabedit using fchk files from gaussian on Windows XP

The goal is to show how to make pretty figures of MOs from Gaussian computations using GabEdit.

This is a write-up for a colleague whom is using Windows. The GabEdit bits work pretty much the same way on linux, although the povray stuff is a lot easier (in my opinion anyway) on that OS.

Anyway, Windows it is.

---

A. Optional: Install Gaussian (G09W)
If you have G09W you can install it. To use the formchk utility, add the location of your gaussian installation to path.

Do this by right-clicking on My Computer, select the Advanced tab, click on Environment Variables, and then click on the Path entry under System Variables. Click Edit, and add

;c:\g09w

if that's the correct path for you. The ";" is just to separate your entry from the previous one. If you need to add another entry afterwards (e.g. for povray), put a ";" after c:\g09w



Either way it's optional, as you will likely have access to formchk on the same computer/system that you ran your job on and can do the chk -> fchk conversion there.

Also, using formchk on a windows computer on a .chk file from a linux computer might not work (which was my case)

An alternative is to set up a smallish linux VM: https://verahill.blogspot.com.au/2016/08/631-small-linux-installation-in.html

---

B. Optional: Add option to open terminal using mouse
Following this: https://stackoverflow.com/questions/378319/windows-explorer-command-prompt-here/379804

Create a plain text file called context.reg



Edit it and put the following in it:

Windows Registry Editor Version 5.00

[HKEY_LOCAL_MACHINE\SOFTWARE\Classes\Folder\shell\cmd]
@="Open Command Prompt Here"
[HKEY_LOCAL_MACHINE\SOFTWARE\Classes\Folder\shell\cmd\command]@="cmd.exe /k pushd %L"

Like so...

Run it (i.e. double-click):

Yes!

You will now have an "Open Command Prompt Here" option when right-clicking on a folder in the file explorer.

Open terminal here

Why do this? Because it makes it easier to run either formchk or povray in the directory where you have your files.

---

C. Install gabedit
Got to https://sites.google.com/site/allouchear/Home/gabedit/download and get "The Gabedit installer for Windows (OpenGL)"

Run. It's that easy.

When running gabedit for the first time you may be asked lots of questions. You can probably just accept all the defaults. We'll change some of them below.

---

D. Install povray
Download http://www.povray.org/redirect/www.povray.org/ftp/pub/povray/Official/povwin-3.7-agpl3-setup.exe and also accept to install the editor dlls at the end of the installation.

Download povray

Once that's done, add povray to path. Go to My Computer, right-click, click on the Advanced tab, select Environment Variables, click on path under System Variables and then Edit. Add

;C:\Program Files\POV-Ray\v3.7\bin

Like So

---

---

Let's get going!

---

1. Optional: Convert chk to fchk using formchk
Didn't work for me, but you can always try (I convert the files under linux)

Either way, you must have the fchk file before going to gabedit

---

2. Open fchk file in gabedit
You need to do a bit of setting up:

Set path and command for povray






To look at the fchk file:
a) Go to Display

Display Geometry/Orbitals/Density/Vibration

b) unset dipole (because otherwise you'll have an annoying arrow in your figure)

Untick dipole

c) Set desired atom colours.
In particular, I change N from dark blue to light blue so that it doesn't interfere with the orbital colours.

Atom colour

d) Set backgrounds (povray and screen) to white

Povray background

Window background

e) load fchk file

Go to Orbitals/Gaussian fchk

Open your fchk file

f) Set ball/stick parameters

Set Ball and Stick parameters

Apply

g) select level/MO

Note that you get the alpha OR beta orbitals. For a spin restricted system these are the same (apart from signage)

You should edit number of points. 65 is the default. 80 is slow but manageable. 90 is pushing it. Higher = smoother.

Computing...

"Get Isovalue" almost never works properly. Try 0.04-0.06.

Et voilá!

To select a different MO, go back to orbitals and click on Selection.

Selection

h) Edit surfaces
You can change colours and transparency

Transparency is an option

Adjust the transparency/opacity here

i) Export pov
To make a publication-worthy figure, export as pov

Export as Povray

Clicking 'Run' occasionally leads to crash, so it's safer to click on save. Make sure that the background in white.

Edit the command to +A0.01 instead of +A0.3 for less 'pixelated' figure on zooming in. Takes a bit longer to render though. The H/W values depend on your window shape.

Select save, not Run

---

3. Render povray
You now have a pov file (example.pov). The most straight-forward way to render it is to open the command prompt (terminal) in the same folder and run

pvengine +A0.01 +H604 +W620 example.pov

Rendering

assuming that the H/W values are the ones listed (or multiples of them) when you saved the pov file (it'll be listed in example.bat in the same folder).

And it looks like this:

A basic POV-Ray generated png

You can also run the .bat file associated with the pov file, but note that you'll have to put paths in "":

You need ""s

You can also render in the povray editor (open the povray with the editor):

Povray editor

Render Settings

The issue becomes the default H/W (resolution) values. Override using command line options: +HXXX +WYYY

You may have to associate .pov files with pvengine manually:

Edit -- but you'll be told that there's no program

Select a program to open

pvengine

Always use the selected program

337. Modifying Nwchem 6.1.1 to work with GabEdit

Karol Strutynski left the following comment on a post about NWChem and Gabedit:

Hello,
I have one important comment:
The vectors coefficients in the nwchem output are incomplete!
The default behaviour of nwchem is to print 10 first coefficients with value bigger than 0.15. For systems with many atoms it is not enough, usually its not even close.

This behaviour is hard-coded in the nwchem source.
To change this you must search each instance of movecs_print_anal in the source code and replace 0.15d0 for smaller value in appropriate calls.
Furthermore you must change one loop in the src/ddscf/movecs_pr_anal.F file and around 200 line there will be loop:
do klo = 0, min(n-1,9), 2
You must increase the range of this loop, for something more reasonable like:
do klo = 0, min(n-1,199), 2

After recompiling the nwchem will print more coefficients and the gabedit will produce more reliable orbitals.

Best regards,
Karol Strutynski

So let's modify NWChem. I'll be modifying the 27th of June release of NWChem 6.1.1, which you'll obtain as Nwchem-6.1.1-src.2012-06-27.tar.gz from http://www.nwchem-sw.org/index.php/Download.

---

Change the number in red to something smaller (I tried 0.01d0) in the following files:
 /src/ddscf/uhf.F

 146  9611    continue
 147          call movecs_print_anal(basis, ilo, ihi, 0.15d0, g_movecs,
 148      $        'UHF Final Alpha Molecular Orbital Analysis',
 149      $        .true., dbl_mb(k_eval), oadapt, int_mb(k_irs),
 150      $        .true., dbl_mb(k_occ))
 151          call movecs_print_anal(basis, ilo, ihi, 0.15d0, g_movecs(2),
 152      $        'UHF Final Beta Molecular Orbital Analysis',
 153      $        .true., dbl_mb(k_eval+nbf), oadapt, int_mb(k_irs+nmo),
 154      $        .true., dbl_mb(k_occ+nbf)

/src/ddscf/scf_vec_guess.F

506          if (scftype.eq.'RHF' .or. scftype.eq.'ROHF') then
507             call movecs_print_anal(basis, 1,
508      &           nprint, 0.15d0, g_movecs,
509      &           'ROHF Initial Molecular Orbital Analysis',
510      &           .true., dbl_mb(k_eval), oadapt, int_mb(k_irs),
511      &           .true., dbl_mb(k_occ))
512          else
513             nprint = min(nalpha+20,nmo)
514             call movecs_print_anal(basis, max(1,nbeta-20),
515      &           nprint, 0.15d0, g_movecs,
516      &           'UHF Initial Alpha Molecular Orbital Analysis',
517      &           .true., dbl_mb(k_eval), oadapt, int_mb(k_irs),
518      &           .true., dbl_mb(k_occ))
519             call movecs_print_anal(basis, max(1,nbeta-20),
520      &           nprint, 0.15d0, g_movecs(2),
521      &           'UHF Initial Beta Molecular Orbital Analysis',
522      &           .true., dbl_mb(k_eval+nbf), oadapt, int_mb(k_irs+nmo),
523      &           .true., dbl_mb(k_occ+nbf))

/src/ddscf/rohf.F

155          endif
156          call movecs_print_anal(basis, ilo, ihi, 0.15d0, g_movecs,
157      $        'ROHF Final Molecular Orbital Analysis',
158      $        .true., dbl_mb(k_eval), oadapt, int_mb(k_irs),
159      $        .true., dbl_mb(k_occ))

/src/mcscf/mcscf.F

680       if (util_print('final vectors analysis', print_default))
681      $     call movecs_print_anal(basis,
682      $     max(1,nclosed-10), min(nbf,nclosed+nact+10),
683      $     0.15d0, g_movecs, 'Analysis of MCSCF natural orbitals',
684      $     .true., dbl_mb(k_evals), .true., int_mb(k_sym),
685      $     .true., dbl_mb(k_occ))
686 c

/src/nwdft/scf_dft_cg/dft_cg_solve.F

166           call movecs_fix_phase(g_movecs(ispin))
167           call movecs_print_anal(basis, ilo, ihi, 0.15d0,
168      &         g_movecs(ispin),blob,
169      &         .true., dbl_mb(k_eval+(ispin-1)*nbf),
170      &         oadapt, int_mb(k_irs+(ispin-1)*nbf),
171      &         .true., dbl_mb(k_occ+(ispin-1)*nbf))
172         enddo

/src/nwdft/scf_dft/dft_scf.F

1736             call movecs_print_anal(ao_bas_han, ilo, ihi, 0.15d0,
1737      &           g_movecs(ispin),
1738      &           blob,
1739      &           .true., dbl_mb(k_eval(ispin)), oadapt,
1740      &           int_mb(k_ir+(ispin-1)*nbf_ao),
1741      &           .true., dbl_mb(k_occ+(ispin-1)*nbf_ao))

/src/nwdft/scf_dft/dft_mxspin_ovlp.F

186       call movecs_print_anal(basis,int_mb(k_non),int_mb(k_non)
187      & ,0.15d0,g_alpha,'Alpha Orbitals without Beta Partners',
188      &   .false., 0.0 ,.false., 0 , .false., 0 )
189 c
190       if (nct.GE.2) then
191       do i = 2,nct
192       ind = int_mb(k_non+i-1)
193       call movecs_print_anal(basis,ind,ind
194      & ,0.15d0,g_alpha,' ',
195      &   .false., 0.0 ,.false., 0 , .false., 0 )
196       enddo
197       endif

352 c
353        call movecs_print_anal(basis, 1, nalp, 0.15d0, g_ualpha,
354      & 'Alpha Orb. w/o Beta Partners (after maxim. alpha/beta overlap)',
355      &   .false., 0.0 ,.false., 0 , .false., 0 )
356 c

Otherwise once could presumably edit the header in ./src/ddscf/movecs_pr_anal.F directly and substitute thresh. At a minimum you should edit that file according to Karol's instructions: change the number in red below to e.g. 199.

/src/ddscf/movecs_pr_anal.F

198             do klo = 0, min(n-1,9), 2
199                khi = min(klo+1,n-1)
200                write(LuOut,2) (
201      $              int_mb(k_list+k)+1,
202      $              dbl_mb(k_vecs+int_mb(k_list+k)),
203      $              (byte_mb(k_tags+int_mb(k_list+k)*16+m),m=0,15),
204      $              k = klo,khi)
205  2             format(1x,2(i5,2x,f12.6,2x,16a1,4x))
206             enddo

---

Compilation
At this point you should be able to follow post 242. Briefly: Compiling NWChem 6.1.1 with Python on Debian Testing (Wheezy) and compile nwchem with python etc. Don't forget to edit /src/config/makefile.h for python support as shown in that post. Once you're done with that you can compare the GabEdit plots with and without the modification.

Alternatively, if you're simply making changes to a copy of nwchem that you've compiled before, you can speed thing up by a factor of ca 300 by following this post:
http://verahill.blogspot.com.au/2013/04/380-modifying-nwchem-code-without-full.html

---

The difference:
I ran a job on benzene as described in post 281. Visualising NWChem output with GabEdit. I chose to run use the ELF (electron localisation function) on output from the unmodified and modified nwchem binaries. It's a pretty big difference:

Original

Modified

281. Visualising NWChem output with GabEdit

Update: please read Karol's comment below. I will put a link here once I've written up a post on how to modify nwchem.

Update 2: Here's the post: http://verahill.blogspot.com.au/2013/02/3xx-modifying-nwchem-611-to-work-with.html .The conclusion is that you MUST edit nwchem. Luckily, it's easy.

Original post:
I've never liked Gabedit much (looks a bit dated, tries to do 'too much') -- until today. Suddenly I have a newfound respect for the developer(s) behind it. It actually doesn't try to do 'too much' -- it simply does A LOT, and actually does it in a pretty transparent way.

Long story short -- you can do things with gabedit which you can't do (easily) with ECCE, and as such it has become an important ally. Besides, it's always nice to have alternatives.

GabEdit is in the Debian repos.

Running your calculations
There are some restrictions"
1. NOTE: you must run your nwchem job with explicit basis sets (i.e. entered as text) -- to do that in ECCE tick the box as shown in the figure below. If you're running 'pure' nwchem, you (probably) have to cut and paste from the basis set directory -- see e.g. section 7.2 here. It's a minor convenience for gaining access to what GabEdit has to offer.

2. You can only open Single point/Energy calculations i.e. Optimizations won't work. So do a single point calculation on your optimized structure.

3. Also, you need to rename/copy your output file so that it ends with .out.

gabedit won't read it otherwise

GabEdit
It's fairly straightforward -- just point and click. One thing which you will want to play with are the iso-surface settings. The defaults are rarely good.

Anyway, I'll let the screenshots do the talking:

Go straight to the Output viewer -- Geometry/Orbital/Density

Click on the M, or right-click anywhere in the window, and load your renamed nwchem output file.

Here's triplet oxygen. The alpha, beta orbitals are listed in the right window

You can do electron localisation

Look at spin density (the unpaired electrons are in the anti-bonding  pi orbitals)

Contour plots are neat -- here showing spin density

Electrostatic potential. 

There's a lot to explore. GabEdit can obviously also prepare and submit jobs, but I'm happy with ECCE in this respect, and content with using GabEdit for post-processing.