21 March 2017

634. Compiling EQ3/6 on Debian

EQ3/6, a program for calculation of geochemical equilibrium in aqueous systems and which can be downloaded here, has a spectacularly unhelpful Install.txt file.

Sadly, we do not provide specific instructions on how to compile and link
the software to build the executables. You have to figure that out yourself
for your specific platform. It will require some number of steps. Also, we 
do not provide makefiles. It is not necessary to use make files, although
some people like them. If all this seems to exceed your knowledge and capacity,
seek the assistance of a qualified Linux or Unix maven who knows how to do
these things on your platform. 


Here I describe how I went about to compile it.


1. Download EQ36_80a_linux.zip from https://www-gs.llnl.gov/about/energy-security/geochemistry


2. Extract it. Move the EQ3_6v8.0a folder to ~
You can put it wherever you want, but then need to edit eq36cfg, and you'll have to deal with csh not handling long paths that well. It's just easier to give up and put EQ3_6v8.0a under $HOME


3. Time to extract further, compile and set up
cd ~/EQ3_6v8.0a
ls *.tar |xargs -I {} tar xvf {}
find -name "*.gz"|xargs -I {} gunzip {}

Time to look at what we have. There are .f files in
xcon6/src
eq3nr/src
eqlib/src
eqlibu/src
eq6/src
xcon3/src
eqpt/src
eqlibg/src

A lot of these are subroutines. We'll compile all files as modules and make static libraries
cd eqlibg/src
gfortran -c *.f
ar cr eqlibg.a *.o
cd ../../eqlibu/src
gfortran -c *.f
ar cr eqlibu.a *.o
cd ../../eqlib/src
gfortran -c *.f
ar cr eqlib.a *.o

cd ../../eqpt/src
sh makelinks
gfortran *.f *.a -o eqpt
cp eqpt ../../bin/
cd ../../xcon3/src
sh makelinks
gfortran *.f *.a -o xcon3
cp xcon3 ../../bin/
cd ../../xcon6/src
sh makelinks
gfortran *.f *.a -o xcon6
cp xcon6 ../../bin/
cd ../../eq3nr/src
sh makelinks
gfortran *.f *.a -o eq3nr
cp eq3nr ../../bin/
cd ../../eq6/src
sh makelinks
gfortran *.f *.a -o eq6

You'll get
eq6.f:16.9: use mod6pt 1 Fatal Error: Can't open module file 'mod6pt.mod' for reading at (1): No such file or directory eqcalc.f:77.9: use mod6pt 1 Fatal Error: Can't open module file 'mod6pt.mod' for reading at (1): No such file or directory eqphas.f:129.9: use mod6pt 1 Fatal Error: Can't open module file 'mod6pt.mod' for reading at (1): No such file or directory eqshel.f:146.9: use mod6pt 1 Fatal Error: Can't open module file 'mod6pt.mod' for reading at (1): No such file or directory exivar.f:121.9: use mod6pt 1 Fatal Error: Can't open module file 'mod6pt.mod' for reading at (1): No such file or directory
Don't despair, just run again:

gfortran *.f *.a -o eq6
cp eq6 ../../bin/

You should now have a bin/ that looks like
eq3nr eq6 eqpt readme.txt runeq3 runeq6 runeqpt xcif3 xcif6 xcon3 xcon6
cd ~/EQ3_6v8.0a
cat eq36cfg >> ~/.cshrc
cd scripts/
ln -s runeq36 runeq3
ln -s runeq36 runeq6
sudo mkdir /usr/tmp
sudo chmod ugo+rw /usr/tmp
cd ../db/
csh ../scripts/runeqpt all
------------------------------------ All done

4. Test
cd ~/EQ3_6v8.0a/3tlib_fmt
csh ../scripts/runeq36 fmt ../db/data0 c4pgwbN2.3i
[..] Done. Optimization ended outside requested limits. Starting hybrid Newton-Raphson iteration. Done. Hybrid Newton-Raphson iteration converged in 56 iterations. * Warning - (EQ3NR/eq3nr) The calculated TDS of 3.44395E+05 mg/L differs from the input file value of 3.54016E+05 mg/L by more than 1%. The calculated value will be used in subsequent calculations. The pickup file has been written. No further input found. Start time = 13:57:15 21Mar2017 End time = 13:57:15 21Mar2017 Run time = 0.586E-01 seconds Normal exit Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG ------------------------------------ The following output files were written: c4pgwbN2.3o c4pgwbN2.3p ------------------------------------ The following input files were run without generating any EQ3NR error messages: c4pgwbN2.3i ------------------------------------ All done
The question here is whether the error is indicative of a serious issue with the compiled software.

I reran as
csh ../scripts/runeq3 fmt c4pgwbN2.3i |tee c4pgwbN2.3o

and compared
diff c4pgwbN2.3o Outputs/c4pgwbN2.3o
18c18 < Run 14:12:44 21Mar2017 --- > Run 13:18:06 07Sep2011 870,871c870,871 < -0.0000 per cent of the total charge < -0.0000 per cent of the mean charge --- > 0.0000 per cent of the total charge > 0.0000 per cent of the mean charge 1235,1236c1235,1236 < Start time = 14:12:44 21Mar2017 < End time = 14:12:44 21Mar2017 --- > Start time = 13:18:06 07Sep2011 > End time = 13:18:06 07Sep2011 1238c1238 < Run time = 0.586E-01 seconds --- > Run time = 0.781E-01 seconds

The results are the same.We're good.

Everything put together as a script:
mkdir ~/tmp cd ~/tmp wget http://www-gs.llnl.gov/content/assets/docs/energy/EQ36_80a_Linux.zip unzip EQ36_80a_Linux.zip cd Linux mv EQ3_6v8.0a ~/ rm -rf Linux cd ~/EQ3_6v8.0a ls *.tar |xargs -I {} tar xvf {} find -name "*.gz"|xargs -I {} gunzip {} cd eqlibg/src gfortran -c *.f ar cr eqlibg.a *.o cd ../../eqlibu/src gfortran -c *.f ar cr eqlibu.a *.o cd ../../eqlib/src gfortran -c *.f ar cr eqlib.a *.o cd ../../eqpt/src sh makelinks gfortran *.f *.a -o eqpt cp eqpt ../../bin/ cd ../../xcon3/src sh makelinks gfortran *.f *.a -o xcon3 cp xcon3 ../../bin/ cd ../../xcon6/src sh makelinks gfortran *.f *.a -o xcon6 cp xcon6 ../../bin/ cd ../../eq3nr/src sh makelinks gfortran *.f *.a -o eq3nr cp eq3nr ../../bin/ cd ../../eq6/src sh makelinks gfortran *.f *.a -o eq6 gfortran *.f *.a -o eq6 cp eq6 ../../bin/ cd ~/EQ3_6v8.0a cat eq36cfg >> ~/.cshrc cd scripts/ ln -s runeq36 runeq3 ln -s runeq36 runeq6 sudo mkdir /usr/tmp sudo chmod ugo+rw /usr/tmp cd ../db/ csh ../scripts/runeqpt all

31 January 2017

633. Multiwfn on debian, with gaussian

Not much to say. I wanted to do some Bader-type AIM analyses, and I don't have AIM2000/AIMPAC. Multiwfn can do this, and a lot more.


Setting up multiwfn
Download the precompiled Multiwfn binaries here: https://multiwfn.codeplex.com/

Extract the binary and the copy the directory to .e.g /opt

$ ls /opt/Multiwfn_3.3.9_bin_Linux/
examples  Multiwfn  settings.ini

Edit your ~/.bashrc and add:
export KMP_STACKSIZE=64000000
export Multiwfnpath=/opt/Multiwfn_3.3.9_bin_Linux

When I hit 0 to visualize something, the plot screen is empty. Check/uncheck a box in the window and the plot will show.
Generating wfx files in gaussian:
If you've already run a job that has generated a .chk file, then run a simple quick job like this:
%nprocshared=2 %mem=500Mb %chk=test.chk #P chkbasis geom=allcheck guess=(check,only) output=wfx pop=(none) test.wfx

To analyse the .wfx file using multiwfn, look at the examples in the manual: https://www.codeplex.com/Download?ProjectName=multiwfn&DownloadId=1607667
The strength of multiwfn is that it can do a LOT. That's also its weakness, in the sense that the menu options are a bit overwhelming. Following the examples in the manual is probably mandatory.

12 January 2017

632. How to set up a few different types of jobs in Gaussian 09 (D/E)

I haven't posted anything for a long time. Partly because with each post I have one fewer issue to solve in the future, and partly because I'm busy, which is a good thing to be if you're an academic.

Either way, I've been doing a lot of fun computational stuff during the past year. Because I have the memory of a goldfish (not as bad as they claim) I'll post some of the techniques here. Nothing fancy, but still handy stuff.

Let's start with the easiest. I'll also have to disappoint by not actually going through the inputs (I don't have the time) -- but at least it gives you something to dissect.


Paramagnetic systems
You should always do a fragment guess when doing computations on a paramagnetic system. There used to be pretty decent example on the gaussian website, but the link's dead.
So, here's an example using MnII(H2O)72+:
%nprocshared=6
%Mem=800000000
%Chk=Mn_heptaaqua_fragmentguess-1.chk
#P uPBE1PBE/def2svp 5D 7F Punch=(MO) Pop=() guess=(fragment=8,only) nosymm

Mn_heptaaqua_fragmentguess

2 6 2 6 0 1 0 1 0 1 0 1 0 1 0 1 0 1! charge and multiplicity
 Mn(fragment=1)    -0.00562000     -0.159820     -0.00534000
 O(fragment=2)    0.806850     2.92945     0.648110
 H(fragment=2)    0.517800     2.22992     1.19065
 H(fragment=2)    0.710610     2.64504     -0.251480
 O(fragment=3)    -2.27505     2.01043     -0.829720
 H(fragment=3)    -1.61846     1.53635     -1.35764
 H(fragment=3)    -2.16217     1.69798     0.0680700
 O(fragment=4)    1.25719     0.434100     -2.53776
 H(fragment=4)    1.44642     -0.489420     -2.38016
 H(fragment=4)    1.93095     0.912890     -2.06836
 O(fragment=5)    -0.763700     -0.365100     2.64881
 H(fragment=5)    -1.25350     0.434770     2.55458
 H(fragment=5)    0.142800     -0.109660     2.78161
 O(fragment=6)    2.85799     -0.0152600     0.774460
 H(fragment=6)    2.85234     -0.657260     0.0556300
 H(fragment=6)    2.56047     -0.559630     1.50359
 O(fragment=7)    -2.58701     -1.35998     -0.392700
 H(fragment=7)    -2.33975     -1.98479     0.179520
 H(fragment=7)    -2.25366     -1.27355     -1.16844
 O(fragment=8)    0.506670     -2.98431     -0.354720
 H(fragment=8)    0.314610     -2.83412     0.445580
 H(fragment=8)    0.228450     -2.74784     -1.07163

--Link1--
%nprocshared=3
%Mem=800000000
%Chk=Mn_heptaaqua_fragmentguess-1.chk
#P uPBE1PBE/chkbasis 5D 7F Opt=() SCRF=(PCM,Solvent=water)  Punch=(MO) Pop=() guess=read geom=allcheck nosymm
The keys here are: 1) everything is divided into fragments (Mn and each H2O are individual fragments). 2) The spin state for the total system is given, followed by the multiplicities and charges of each fragment and 3) the first calculation is just used to generate a fragment guess ( guess=(fragment=8,only) ). 4) The first job is then immediately followed by a second job defined in link1, and it reads in the fragment guess, the geometry and the basis set from the first job. NOTE: fragment guesses can't use SCRF.

Freezing atoms
This is pretty basic, but useful. The example uses Mg(H2O)62+.
%nprocshared=6
%Mem=800000000
%Chk=Mg6aq_b3pw91.chk
#P rPBE1PBE/gen 5D Opt=(readopt)  Punch=(MO) Pop=() 

Mg6aq_b3pw91

2 1 ! charge and multiplicity
 Mg     0.00000     0.00000     0.00000
 O     0.00000     2.09000     0.00000
 H     -0.691981     2.65500     -0.691981
 H     0.691981     2.65500     0.691981
 O     1.47785     2.22045e-16     1.47785
 H     1.87737     0.978609     1.87737
 H     1.87737     -0.978609     1.87737
 O     -1.47785     -1.11022e-16     1.47785
 H     -1.18539     7.33956e-09     2.56935
 H     -2.56935     -7.33957e-09     1.18539
 O     0.00000     -2.09000     0.00000
 H     -0.691981     -2.65500     0.691981
 H     0.691981     -2.65500     -0.691981
 O     -1.47785     2.22045e-16     -1.47785
 H     -1.87737     -0.978609     -1.87737
 H     -1.87737     0.978609     -1.87737
 O     1.47785     -1.11022e-16     -1.47785
 H     1.18539     -2.20187e-08     -2.56935
 H     2.56935     2.20187e-08     -1.18539

noatoms atoms=Mg,O

Mg H O 0
6-31G
****

Keys: 1) use opt=(readopt) to make gaussian look for more instructions immediately noatoms makes all atoms inactive. Then make Mg and O active again using 'atoms'. You can also give a numerical range, e.g. atoms=1-2,5,8 etc. You can also use notatoms, which specifies atoms to make inactive (instead of noatoms which makes all inactive).

I defined the basis set at the bottom only to show what order things come in i.e. readopt input comes before basis gen.

NOTE: if you do opt=(readopt) and freq=() in the same job those atoms will be frozen in the Hessian

Mixing basis sets
What if you have a large cluster and want to use a triple-zeta basis set for some oxygens, and a double-zeta for others that are less 'important'?
%nprocshared=6
%Mem=800000000
%Chk=Mg6aq_b3pw91.chk
#P rPBE1PBE/gen 5D Opt=(readopt)  Punch=(MO) Pop=() 

Mg6aq_b3pw91

2 1 ! charge and multiplicity
 Mg     0.00000     0.00000     0.00000
 O     0.00000     2.09000     0.00000
 O     -1.47785     -1.11022e-16     1.47785
 O     0.00000     -2.09000     0.00000
 O     -1.47785     2.22045e-16     -1.47785
 O     1.47785     2.22045e-16     1.47785
 O     1.47785     -1.11022e-16     -1.47785
 H     -0.691981     2.65500     -0.691981
 H     0.691981     2.65500     0.691981
 H     1.87737     0.978609     1.87737
 H     1.87737     -0.978609     1.87737
 H     -1.18539     7.33956e-09     2.56935
 H     -2.56935     -7.33957e-09     1.18539
 H     -0.691981     -2.65500     0.691981
 H     0.691981     -2.65500     -0.691981
 H     -1.87737     -0.978609     -1.87737
 H     -1.87737     0.978609     -1.87737
 H     1.18539     -2.20187e-08     -2.56935
 H     2.56935     2.20187e-08     -1.18539

Mg 0
6-311G
****
O 0
def2svp
****
8-12 0
6-31G
****
13-19 0
6-311+G(2d,2p)
****



In this case Mg uses 6-311G, O uses def2svp, and some of the H use 6-31G while others use 6-311+G(2d,2p). The input obviously doesn't make any chemical sense, but it works as an example. To add ECP, put in an empty space after the last basis set, and then define the ECPs like you would basis sets (but use e.g. def2 instead of def2tzvp, since def2 contains the ecps).

Partial Charges in DFT input
I'm not 100% sure about this one, but I've checked against nwchem and it gives the same output as using 'charge' in nwchem input.
%nprocshared=6
%Mem=800000000
%Chk=Mg6aq_b3pw91.chk
#P rPBE1PBE/GEN 5D Opt=()  Punch=(MO) Pop=() 

Mg6aq_b3pw91

2 1 ! charge and multiplicity
 Mg    0.00000     0.00000     0.00000
 O     0.00000     2.09000     0.00000
 O     1.47785     2.22045e-16     1.47785
 O     -1.47785     -1.11022e-16     1.47785
 O     0.00000     -2.09000     0.00000
 O     -1.47785     2.22045e-16     -1.47785
 O     1.47785     -1.11022e-16     -1.47785
 H(znuc=.91666666666666666666)   -0.691981     2.65500     -0.691981
 H(znuc=.91666666666666666666)   0.691981     2.65500     0.691981
 H(znuc=.91666666666666666666)   1.87737     0.978609     1.87737
 H(znuc=.91666666666666666666)   1.87737     -0.978609     1.87737
 H(znuc=.91666666666666666666)   -1.18539     7.33956e-09     2.56935
 H(znuc=.91666666666666666666)   -2.56935     -7.33957e-09     1.18539
 H(znuc=.91666666666666666666)   -0.691981     -2.65500     0.691981
 H(znuc=.91666666666666666666)   0.691981     -2.65500     -0.691981
 H(znuc=.91666666666666666666)   -1.87737     -0.978609     -1.87737
 H(znuc=.91666666666666666666)   -1.87737     0.978609     -1.87737
 H(znuc=.91666666666666666666)   1.18539     -2.20187e-08     -2.56935
 H(znuc=.91666666666666666666)   2.56935     2.20187e-08     -1.18539

Mg O H 0
6-31G
****

Again, it's an artificial example. I've set the charges of the H to 11/12ths, so that the charge of the entire complex is 1, and not 2+ anymore. Confusingly, I can only get it to work properly if I set the charge to 2+ in the input i.e. ignoring the changed charges for different atoms. Either way, it gives the same result as nwchem (although in nwchem I'd set the total charge to +1 in recognition of the change atomic charges).