Showing posts with label computational. Show all posts
Showing posts with label computational. Show all posts

24 July 2015

614. SIESTA with MPI and acml on debian jessie

One of my students might be using SIESTA for some simulations, and a first step towards that is to set it up on my cluster.

This isn't an optimised build -- right now I'm just looking at having a simple parallell build that runs.

I had a look at and
NOTE: don't use the int64 acml or openblas BLAS libs, or you'll get SIGSEV due to invalid memory reference when running. NWChem is the complete opposite, and for some reason both the int64 and regulat acml libs have the same names. Not sure how that's supposed to work out on a system with nwchem, which needs the int64 libs.

See here for acml on debian. I've got /opt/acml/acml5.3.1/gfortran64_int64/lib in my /etc/ on behalf of nwchem.
 Being lazy, I opted for the debian scalapack and libblacs packages:
sudo apt-get install libscalapack-mpi-dev libblacs-mpi-dev libopenmpi-dev

To get the link to the SIESTA code, go to

Then, if you're an academic, you can do:
sudo mkdir /opt/siesta
sudo chown $USER /opt/siesta
cd /opt/siesta
tar xvf siesta-3.2-pl-5.tgz
cd siesta-3.2-pl-5/Obj
sh ../Src/
*** Compilation setup done. *** Remember to copy an arch.make file or run configure as: ../Src/configure [configure_options]
../Src/./configure --help
`configure' configures siesta 2.0 to adapt to many kinds of systems. Usage: ./configure [OPTION]... [VAR=VALUE]... [..] Installation directories: --prefix=PREFIX install architecture-independent files in PREFIX [/usr/local] --exec-prefix=EPREFIX install architecture-dependent files in EPREFIX [PREFIX] By default, `make install' will install all the files in `/usr/local/bin', `/usr/local/lib' etc. You can specify an installation prefix other than `/usr/local' using `--prefix', for instance `--prefix=$HOME'. [..] --enable-mpi Compile the parallel version of SIESTA --enable-debug Compile with debugging support --enable-fast Compile with best known optimization flags Optional Packages: --with-PACKAGE[=ARG] use PACKAGE [ARG=yes] --without-PACKAGE do not use PACKAGE (same as --with-PACKAGE=no) --with-netcdf=<lib> use NetCDF library --with-siesta-blas use BLAS library packaged with SIESTA --with-blas=<lib> use BLAS library --with-siesta-lapack use LAPACK library packaged with SIESTA --with-lapack=<lib> use LAPACK library --with-blacs=<lib> use BLACS library --with-scalapack=<lib> use ScaLAPACK library [..]
../Src/./configure --enable-mpi
checking build system type... x86_64-unknown-linux-gnu checking host system type... x86_64-unknown-linux-gnu [..] checking for mpifc... no checking for mpxlf... no checking for mpif90... mpif90 checking for MPI_Init... no checking for MPI_Init in -lmpi... yes [..] checking for sgemm in /opt/openblas/lib/ yes checking LAPACK already linked... yes checking LAPACK includes divide-and-conquer routines... yes configure: using DC_LAPACK routines packaged with SIESTA due to bug in library. Linker flag might be needed to avoid duplicate symbols configure: creating ./config.status config.status: creating arch.make
Edit arch.make:
# # This file is part of the SIESTA package. # # Copyright (c) Fundacion General Universidad Autonoma de Madrid: # E.Artacho, J.Gale, A.Garcia, J.Junquera, P.Ordejon, D.Sanchez-Portal # and J.M.Soler, 1996- . # # Use of this software constitutes agreement with the full conditions # given in the SIESTA license, as signed by all legitimate users. # .SUFFIXES: .SUFFIXES: .f .F .o .a .f90 .F90 SIESTA_ARCH=x86_64-unknown-linux-gnu--unknown FPP= FPP_OUTPUT= FC=mpif90 RANLIB=ranlib SYS=nag SP_KIND=4 DP_KIND=8 KINDS=$(SP_KIND) $(DP_KIND) FFLAGS=-g -O2 FPPFLAGS= -DMPI -DFC_HAVE_FLUSH -DFC_HAVE_ABORT LDFLAGS= ARFLAGS_EXTRA= FCFLAGS_fixed_f= FCFLAGS_free_f90= FPPFLAGS_fixed_F= FPPFLAGS_free_F90= BLAS_LIBS=-L/opt/acml/acml5.3.1/gfortran64/lib -lacml LAPACK_LIBS= BLACS_LIBS=-L/usr/lib -lblacs-openmpi -lblacsCinit-openmpi SCALAPACK_LIBS=-L/usr/lib -lscalapack-openmpi COMP_LIBS=dc_lapack.a NETCDF_LIBS= NETCDF_INTERFACE= MPI_LIBS= -L/usr/lib/openmpi/lib -lmpi -lmpi_f90 LIBS=$(SCALAPACK_LIBS) $(BLACS_LIBS) $(LAPACK_LIBS) $(BLAS_LIBS) $(NETCDF_LIBS) $(MPI_LIBS) -lpthread #SIESTA needs an F90 interface to MPI #This will give you SIESTA's own implementation #If your compiler vendor offers an alternative, you may change #to it here. MPI_INTERFACE=libmpi_f90.a MPI_INCLUDE=. #Dependency rules are created by autoconf according to whether #discrete preprocessing is necessary or not. .F.o: $(FC) -c $(FFLAGS) $(INCFLAGS) $(FPPFLAGS) $(FPPFLAGS_fixed_F) $< .F90.o: $(FC) -c $(FFLAGS) $(INCFLAGS) $(FPPFLAGS) $(FPPFLAGS_free_F90) $< .f.o: $(FC) -c $(FFLAGS) $(INCFLAGS) $(FCFLAGS_fixed_f) $< .f90.o: $(FC) -c $(FFLAGS) $(INCFLAGS) $(FCFLAGS_free_f90) $<
cd ../
ln -s Obj/siesta siesta

I added /opt/siesta/siesta-3.2-pl-5 to $PATH.

To test, edit /opt/siesta/siesta-3.2-pl-5/
6 #SIESTA=../../../siesta 7 SIESTA=mpirun -n 2 ../../../siesta
cd /opt/siesta/siesta-3.2-pl-5/Tests/h3po4_2
export LD_LIBRARY_CONFIG=/opt/acml/acml5.3.1/gfortran64/lib 
>>>> Running h3po4_2 test... ==> Copying pseudopotential file for H... ==> Copying pseudopotential file for O... ==> Copying pseudopotential file for P... ==> Running SIESTA as mpirun -n 2 ../../../siesta ===> SIESTA finished successfully

Also, look at work/h3po4_2.out:
* Running on    2 nodes in parallel
>> Start of run:  24-JUL-2015  21:58:13

                           *  WELCOME TO SIESTA  *       

reinit: Reading from standard input
elaps:  optical           1       0.000       0.000     0.00
>> End of run:  24-JUL-2015  21:58:20

05 July 2012

202. Reproducing the paper: computing organic reduction potentials

This post is not altogether finished yet -- will be updated with additional results as they come in

This isn't as much as a straight reproduction as a test of the authors' thesis that the SCFE+solvation energies can be used vs the results from using a classic adiabatic thermodynamic cycle.

The paper by Speelman and Gillmore can be found here:

Again, I'm doing this in large part to deepen my understanding of the practical aspects of computational chemistry. It's not a critical evaluation of the authors' approach -- I don't possess that kind of expertise.

0. Exploring solvation
Because NWChem implements one implicit solvation model (COSMO) and Gaussian implements three, it may be worth quickly exploring what they yield in terms of predicting solvation energies. In particular for processes like reduction or oxidation, solvation is highly important.

The structure of benzene was optimised at 3-21G in vacuo in both software packages. 3-21G was used throughout since it's a matter of making a crude comparison. The solvent was set to be water for all calcs. Note that benzene is neutral, so the solvation energy will be much smaller than for e.g. a cation.

Programme   Solvation model Energy (hartree)    Relative energy to gas phase
NWChem     none                      -230.9757626747    0
NWChem     COSMO               -230.9820692591    3.96 kcal/mol
G09               none                     -230.97576997        0
G09               PCM                    -230.9795649          2.38 kcal/mol
G09               cPCM                  -230.9796073          2.41 kcal/mol
G09               iPCM                   -230.9810728          3.33 kcal/mol

Basically, the gas phase energies are very similar for both G09 and NWChem. The differences lie in the solvation model energies, which iPCM and COSMO yielding the values closest to each other. Ultimately, the energy of solvation spans a fairly small region.

What about frequency calculations?
Well, without PCM I get an enthalpy correction of 0.106934 and with PCM I get 0.106976, a difference of .026355126 hartree. We're talking about a precision beyond that of the method itself. The entropy varies similarly: 68.754 w/o PCM, 68.764 cal/MolK with or ca 0.003 kcal/mol at 298.15 K.

1. Thermodynamic cycle of compound 5 (1,2,4,5-tetracyano-benzene)

A. G09 was used to optimise the structure first using b3lyp/3-21G, then using b3lyp/6-31+G* while applying a solvation model using PCM (water). The optimised structure was then used with IPCM. Also, the frequencies were calculated in vacuo.

(species:   gas phase E/IPCM E/Free Energy corr)
Neutral:   -601.2211700670/-601.2223222/0.057097
Reduced: -601.3655493120/-601.3680354/0.054093

Absolute potential: 
[(-601.3680354+0.054093)-(-601.2223222+0.057097)]*27.2107 eV/hartree= -4.0469 eV
The paper we're following gives 4.12 as the SCE absolute potential.
4.0469-4.12=-0.0731 V

Not close at obvious problems is the reference potential though -- is this really applicable to acetonitrile? Our ideal reference potential should be somewhere around  4.05+0.7=4.75 eV in order to give the desired value of -0.7 V.

IUPAC gives SHE in acetonitrile as 4.60 eV (as mentioned here by Davis and Fry). Also, Pavlishchuck and Addison discuss different reference electrodes in acetonitrile in Inorg. Chim. Acta 200, 298, 97-102 and recommend subtracting 244 mV from SHE to yield the SCE potential. This still leaves us with 4.60-0.24=4.36 V. That in turn gives 4.05-4.36=-0.31 V as the potential. An improvement, certainly, but still 360 mV off.

If we use the old trick of optimising in one level, but using single-point and solvation energies from another level of theory (see section 3 for the origin of the values):

cpcm: [(-601.513973+0.054093)-(-601.3684665+0.057097)]*27.2107= -4.0411/(-1 e)
ipcm: [(-601.5166797+0.054093)-(-601.3698362+0.057097)]*27.2107=-4.0775/(-1 e)

It hardly changes it at all. Also, it quite clear that solvation effectsdominate over thermochemical correction terms here.

B. NWChem was used to optimise the structure at b3lyp/6-31+G* (gas phase) starting with the b3lyp/3-21G structure obtained from G09. 

(species:   gas phase opt/COSMO E/Enthalpy corr/Entropy Corr)
Neutral:    -601.1996294139/-601.2262889216/66.800 kcal/105.187 cal (times/2 cores: 15 min + 4 min +2h 10 min)
Reduced:  -601.2988217652/-601.373759558/65.498 kcal/105.799 cal (times/2 cores: 15 min + 4 min + 2h 30 min)

Absolute potential:
((-601.373759558+(65.498-298.15*105.799/1000)/627.503)-(-601.2262889216+(66.8-298.15*105.187/1000)/627.503))*27.2107=-4.077 eV

Basically the same as we saw with gaussian (ipcm).

2. The proposed method using 6-311++G**
Instead of using MIDI! I optimised the structures with ub3lyp/6-31+G*/PCM.

(species:  cpcm/ipcm energies in hartree)
Neutral:   -601.3684665/-601.3698362
Reduced: -601.513973/-601.5166797
EQM: -3.959 eV /-3.99571
PotCPCM: (-3.959+4.6067)/(-1.0886)    = -0.5949 eV; reported (comp.) -0.604 V
PotIPCM: (-3.99571+4.6291)/(-1.0508) = -0.603 eV

3. Reality
-0.74 V vs SCE in acetonitrile.

4. Summary
From the computational chemist's point of view anions are difficult. Metals are difficult. Radicals are difficult. What I didn't appreciate is quite how difficult anions really are. There should be no appreciable difference between the computation of an oxidation potential and a reduction potential since they are each other's conjugate-- the difference in these particular examples is only in the nature of the product when starting with a neutral parent species.
Empirically we've been getting excellent results in the lab calculating the oxidation potential of neutral organic compounds, and awful results for the reduction potential. This can't really be blamed on the reference potentials used either since they should be the same for both processes. Instead the challenge must be more fundamental.