676. Mixed basis set with gaussian when continuing from a chk-file

 I'm currently working with computations on paramagnetic species, so being able to read the geometry from the chk-file is important in order to avoid having to define the multiplicities of all the segments of the molecule. 

I also can't use def2-tzvp for the entire molecule since the .chk files during freq cals become too big for my current storage allocation (>500 gb). The work is exploratory at the moment, so compromises are OK.

I first run a calc with a small basis set:

%nprocshared=28
%Mem=24gb
%Chk=I1_hs.chk
#P uPBE1PBE/def2svp 5D 7F Punch=(MO) Pop=() guess=(fragment=8,only) nosymm

intermediate 1, pbe0/def-svp.

0 9 -2 1 0 1 0 1 -3 1 -1 1 0 1 4 4 2 6 
S(fragment=1)          -7.953777000    -0.003180000    -2.912460000
O(fragment=1)          -3.008379000    -0.432531000    -1.322496000
...

--Link1--
%nprocshared=28
%mem=24gb
%chk=I1_hs.chk
#P uPBE1PBE/chkbasis 5D 7F Opt=() SCRF=(PCM,Solvent=water)  Punch=(MO) Pop=() guess=read geom=allcheck nosymm

I then want to optimise the structure obtained at def2-svp using a mixed basis set where I use def2-tzvp for all the cobalt centres, and def2-svp for the organic bits.

%nprocshared=28
%Mem=24gb
%oldChk=I1_hs.chk
%chk=I1_hs_tzvp.chk
#P uPBE1PBE/genecp 5D 7F Opt=() SCRF=(PCM,Solvent=water)  Punch=(MO) Pop=() guess=read geom=allcheck nosymm

C H N O 0
def2svp
****

Co 0
def2tzvp
****

Co 0
def2

The basis set definitions are in blue, and the ecp bit is in red.

657. More on charges in nwchem and gaussian

A now ten-year old paper introduced the concept of Pauling bond-strength conserving terminations (PBS ) in the use of molecular codes for calculations involving extended crystalline systems ('Quantum-Chemical Calculations of Carbon-Isotope Fractionation in CO2(g), Aqueous Carbonate Species, and Carbonate Minerals' by James R. Rustad, Sierra L. Nelmes, Virgil E. Jackson, and David A. Dixon --  see link). The authors used NWChem for the calculations, most likely due to the affiliation between the lead author and PNNL, where NWChem is developed, and where the researchers have been banned from using Gaussian.

I use Gaussian almost exclusively these days, mainly due to how fast it is.

Unfortunately, Gaussian and NWChem behave quite differently when it comes to introduction of specified nuclear charges, so I here compare the two codes in terms of how to set up PBS calculations.

NWChem (6.8):

scratch_dir /scratch
Title "charge"

Start  charge

echo

charge 0

geometry noautosym noautoz units angstrom
 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
 H1     -0.691981     2.65500     -0.691981 charge 0.5 
 H1     0.691981     2.65500     0.691981 charge 0.5 
 H1     1.87737     0.978609     1.87737 charge 0.5 
 H1     1.87737     -0.978609     1.87737 charge 0.5 
 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
end

basis "ao basis" spherical print
  H library "def2-svp"
  Mg library "def2-svp"
  O library "def2-svp"
END

dft
  mult 1
  direct
  XC pbe0
  grid xfine
  mulliken
end

task dft energy   

This gives an energy of -655.860806066326.

Removing the charges for H1 and setting the total charge to +2 gives an energy of -657.044328628867

Gaussian (16.A01):
WRONG:

%nprocshared=6
%Mem=800000000
%Chk=charge.chk
#P GFINPUT rPBE1PBE/def2svp 5D  NoSymm  Punch=(MO) Pop=(full) 

charge

0 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=0.5)     -0.691981     2.65500     -0.691981
 H(znuc=0.5)     0.691981     2.65500     0.691981
 H(znuc=0.5)     1.87737     0.978609     1.87737
 H(znuc=0.5)     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

gives an energy of -655.679686484!

However,

2 1  ! charge and multiplicity

gives an energy of -655.860712881, which is what we want.

Removing the znuc specifications and using

2 1  ! charge and multiplicity

gives an energy of -657.044229333

Keeping the znuc specifications and defining those protons as fragment 2, and the rest of the cluster as fragment 1

 2 1 -2 1 4 1! charge and multiplicity

gives an energy of -655.860712881


Conclusion: 
both NWChem and Gaussian can be made to use PBS, but while you use the intended cluster charge (0) in NWChem, you need to use the unmodified charge (+2) in gaussian.

647. Linda/G16 and ecce + slurm

I had the opportunity to test G16 with Linda on my cluster. Suprisingly, I observed a decent reduction in computational time for geometry optimisation, and a fair reduction for vibrational analysis.

Computational time -- geometry optimisation
1x6 cores: 100%
2x6 cores: 56% (2 nodes)
3x6 cores: 44% (3 nodes)

Computational time -- normal mode computation
1x6 cores: 100%
2x6 cores: 63% (2 nodes)
3x6 cores: 57% (3 nodes)

In no case was the speed up good enough that it actually makes long-term sense to use Linda (need less than 50% for 2 nodes and 33% or less for 3 nodes), but if you're itching for quick results, it might be worth it.

Either way, as part of this I wrote a CONFIG.Linda for ECCE and SLURM:

Gaussian-16: /opt/gaussian/g16a/g16/g16
perlPath: /usr/bin/
qmgrPath: /usr/bin/
xappsPath: /usr/bin/

Slurm {
#!/bin/csh
#SBATCH -p linda
#SBATCH --time=$walltime
#SBATCH --output=slurm.out
#SBATCH --job-name=$submitFile
#SBATCH --nodes=$nodes
#SBATCH --cpus-per-task=$ppn
}

Gaussian-16FilesToRemove{ core }

Gaussian-16Command{
set path = ( /opt/nbo6/bin $path ) 

set nnodes = `scontrol show job $SLURM_JOBID|grep NodeList|sed 's/=/\t/'|gawk '{print $2}'|tail -n 1`
sed -i "2s/^/%LindaWorkers=$nnodes\n/" g16.g16in
sed -i '3s/^/%UseSSH\n/' g16.g16in
sed -i "s/nprocshared.*/nprocshared=$ppn./" g16.g16in

setenv GAUSS_SCRDIR /home/me/scratch
setenv GAUSS_EXEDIR /opt/gaussian/g16a/g16
setenv GAUSS_ARCHDIR /opt/gaussian/g16a/g16/arch
setenv GAUSS_BSDDIR /opt/gaussian/g16a/g16/bsd
setenv GAUSS_LEXEDIR /opt/gaussian/g16a/g16/linda-exe
/opt/gaussian/g16a/g16/g16< $infile > $outfile 
}

The slurm partition linda has all the nodes that has g16.