148. popt, readopt and g09?

Update:
NOTE that readopt WORKS in G09 rev. B and rev. D. It DOESN'T work in G09 rev. A.

readopt is a whole lot easier to use than popt.

Original post:
For various reasons I want to be able to freeze some atoms while including other in my optimisation using gaussion 09. Readopt (http://www.gaussian.com/g_tech/g_ur/k_opt.htm) sounded like a good idea.

The problem
I can't get this to work:

%nprocshared=6
%Chk=reduced_solvate.chk
#P ub3lyp/lanl2dz 6D 10F opt=ReadOptimize Punch=(MO) Pop=() 

reduced_solvate

3 2 ! charge and multiplicity
 Cr     10.0000     10.0000     10.1100
 O     10.0000     10.0000     12.1300
 H     10.0000     10.8000     12.7000
 H     10.0000     9.20000     12.7000
 O     10.0000     7.84000     10.0000
 H     10.0000     7.15000     10.6900
[..]
 H     6.32536     3.63314     12.3506
 H     8.48233     15.9732     5.80567

noatoms atoms=1-19

readopt and rdopt don't work either -- they all give 'syntax error' (C64 flashback!)

 QPErr --- A syntax error was detected in the input line.
 #P ub3lyp/lanl2dz 6D 10F ReadOptimize Pu
                                                  '
 Last state="GCL"

The solution: 

I didn't come up with this: http://ccl.net/cgi-bin/ccl/message-new?2012+04+24+002+raw

Anyway, you can use popt but it takes a bit more preparation of the coordinates.

You can convert an .xyz file quickly by using:

cat molecule.xyz | gawk '{print $1,"-1",$2,$3,$4}'>popt_molecule.xyz

Change 0 to -1 for the atoms you want to optimise.

%nprocshared=6
%Chk=popt_solvate.chk
#Popt ub3lyp/lanl2dz 6D 10F Punch=(MO) scf=(maxcycle=1024)

reduced_solvate popt

3 2 ! charge and multiplicity
Cr 0 10.0000 10.0000 10.1100
O 0 10.0000 10.0000 12.1300
H 0 10.0000 10.8000 12.7000
H 0 10.0000 9.20000 12.7000
O 0 10.0000 7.84000 10.0000
H 0 10.0000 7.15000 10.6900
[..]
H -1 3.62686 7.14993 12.8142
H -1 5.12696 15.5239 9.01700
H -1 6.32536 3.63314 12.3506
H -1 8.48233 15.9732 5.80567

147. Oniom in nwchem -- with a little bit of help from gromacs and openbabel

Example -- I want to do explicit solvent modelling of methanol in water. This is obviously an articifical approach, but generally applicable.

This is a rough approach to doing oniom calculations using nwchem 6.0 -- this is a technical description, not a how-to when it comes to the science.

1. Pre-optimisation
Draw methanol and set up a simple calc using e.g. ecce to pre-optimise the structure with e.g. an implicit solvent model. Here's nwch.nw:

scratch_dir /scratch
Title "pre-oniom"
Start  pre-oniom
echo
charge 0
geometry autosym units angstrom
 C     0.00000     0.00000     0.00000
 H     -0.675500     -0.675500     0.675500
 H     0.675500     -0.675500     -0.675500
 H     -0.675500     0.675500     -0.675500
 O     0.866025     0.866025     0.866025
 H     1.51843     0.213620     1.51843
end
ecce_print /home/me/jobs/jobs/testing/old/pre-oniom/ecce.out
basis "ao basis" spherical print
  H library "6-31+G*"
  O library "6-31+G*"
  C library "6-31+G*"
END
dft
  mult 1
  XC b3lyp
  mulliken
end
driver
  default
end
cosmo
end
task dft optimize
task dft freq

2. Solvation using gromacs

Take the output, nwch.nwout, and use babel to export the optimised structure

babel -inwo nwch.nwout -oxyz molecule.xyz

The next few steps require gromacs:
editconf -f molecule.xyz -o molecule.gro -box 2 2 2
genbox -cp molecule.gro -cs spc216.gro -o solvated.gro
babel -igro solvated.gro -oxyz solvated.xyz

tail -n +3 solvated.xyz > oniom.nw

3. Putting it all together

The only 'trick' is how to define what part of the input belongs to the high level section and what belongs to the low level section -- for a solvation case like this, where whole molecules are treated by one method or another, use model. The last atom to be part of the high level section is the last of the six atoms in methanol, so the keyword is model 6. Atoms 7 to infinity are thus part of the sto-3g part, and atoms 1-6 part of the 6-31g* part.


memory stack 600 mb heap 200 mb global 800 mb

scratch_dir /scratch
Title "oniom"
Start oniom
echo
charge 0

geometry units angstrom
symmetry c1
C         10.12000       10.35000       10.00000
H          9.60000       10.72000       10.89000
H          9.60000       10.72000        9.11000
H         11.14000       10.74000       10.00000

[.]

H          0.46000        2.39000        6.86000

H          1.16000        1.73000        8.18000

end

basis sto-3g spherical

 * library sto-3g

end

basis 6-31g spherical

 * library 6-31g

end

oniom

    model 6

    low  dft basis sto-3g input "dft\; xc\; end"

    high dft basis "6-31g" input "dft\; xc\; end"

end

basis "ao basis" spherical print

  H library "6-31+G*"

  O library "6-31+G*"

  C library "6-31+G*"

END

task oniom

This takes forever to run, but run it does. The memory statement is important -- if the global memory is too small it will crash. Also, be aware that the amount of memory specified is per instance -- if you launch with mpirun -n 3, multiply by 3 to get the amount of memory that needs to be present (here, 2.4 GB physical RAM).

146. Nwchem with openblas

Openblas seems to be the successor of Gotoblas. http://xianyi.github.com/OpenBLAS/
I willingly admit that I am no expert on this blas library stuff, so I may well be doing something obviously wrong.

1. Compiling and installing openblas
sudo mkdir /opt/openblas
sudo chown ${USER} /opt/openblas
mkdir ~/tmp
cd ~/tmp

download from here http://github.com/xianyi/OpenBLAS/tarball/v0.1.1

tar xvf xianyi-OpenBLAS-v0.1.1-0-g5b7f443.tar.gz

 make all BINARY=64 CC=/usr/bin/gcc FC=/usr/bin/gfortran USE_THREAD=0 INTERFACE64=1 1> make.log 2>make.err

make PREFIX=/opt/openblas install
cp lib*.*  /opt/openblas/lib

cd /opt/openblas/lib
rm libopenblas.so.0
ln -s libopenblas.so libopenblas.so.0

do ls /opt/openblas/lib and note what cpu specific file you have such as libopenblas_barcelona-r0.1.1.so or libopenblas_nehalem-r0.1.1.a

edit /etc/ld.so.conf and add

/opt/openblas/lib

2. Compiling nwchem

sudo apt-get install build-essential gfortran python2.7-dev
mkdir /opt/nwchem/

sudo chown ${USER} /opt/nwchem

cd /opt/nwchem

wget http://www.nwchem-sw.org/images/Nwchem-6.0.tar.gz
tar -xvf Nwchem-6.0.tar.gz
cd nwchem-6.0/

For python support, edit src/config/makefile.h and add -lz -lssl to line 1962 (see here: http://verahill.blogspot.com.au/2012/04/adding-python-support-to-nwchem-under.html)

Next, continue
export LARGE_FILES=TRUE
export TCGRSH=/usr/bin/ssh
export NWCHEM_TOP=`pwd`
export NWCHEM_TARGET=LINUX64
export NWCHEM_MODULES="all python"
export PYTHONHOME=/usr
export PYTHONVERSION=2.7
export USE_MPI=y
export USE_MPIF=y
export MPI_LOC=/usr/lib/openmpi/lib
export MPI_INCLUDE=/usr/lib/openmpi/include
export BLASOPT="-L/opt/openblas/lib -lopenblas -lopenblas_barcelona-r0.1.1"
export LIBRARY_PATH=$LIBRARY_PATH:/usr/lib/openmpi/lib
export LIBMPI="-lmpi -lopen-rte -lopen-pal -ldl -lmpi_f77 -lpthread"
cd $NWCHEM_TOP/src
make clean
make nwchem_config
make FC=gfortran

Now edit your ~/.bashrc and add

export PATH=$PATH:/opt/nwchem/nwchem-6.0/bin/LINUX64
export LD_LIBRARY_PATH=/opt/openblas/lib:$LD_LIBRARY_PATH 

and you should be done.

You might find that ecce still won't properly execute nwchem jobs (you don't get an error -- it just doesn't work) -- edit your /etc/ld.so.conf and add a line saying

/opt/openblas/lib

Then do sudo ldconfig and it will work

Time to look into benchmarking...

---

So far I've got the following benchmarks (in seconds) for nwchem when optimising benzene with 6-311++g** using mpirun -n 3.

without -lpthread
_____Internal blas_____openblas_____ATLAS_____acml_____
Ta     249.8                    277.7                 335.9           ------
Be     408.2                    421.7                 568.4        
B      546.4                     579.9                 760.2        

Noticing a mention regarding -lpthread online I compiled with that as part of the mpilibs, and lo and behold:

with -lpthread

_____Internal blas_____openblas_____ATLAS_____acml_____
Ta      201.7                    198.9                233.9           ------
Be      371.1                    314.4                390.3        
B        496.6                    416.8                 521.9

* All run using mpirun -n 3. Be=phenom II X6. B= Athlon II X3. Ta= i7-1600 X4. Intel MKL cost money, acml are free.

This doesn't mean that one library is superior to the other -- it just means that under the conditions I employed and in the presence of any mistake I may have made, these are my observations. Your mileage will vary.