Uncategorized

Dispersion correction for Schrödinger Jaguar Results

I am sharing here the perl scripts that I have used for the calculation of the dispersion corrected energy of Jaguar optimised structures. I confess that I haven’t looked at these scripts for a while and there may (will) be errors. The scripts are certainly very rough and I really ought to tidy them up. The general usage for the script is to call multiple results files through a bash script (below), which then calls the Perl script on each file. The results are then printed to an output file.

At some point, I had intended to make this accept multiple input formats and to act interactively, but as Stefan Grimme has provided his own program, DFTD3, which you can download from here, I am not sure that it is necessary. The main reason for this program is that Stefan’s program accepts XYZ coordinates, and we had multiple files in jaguar format. Also, we wanted to tie a script into Gaussian ’03 in order to carry out dispersion corrected optimisation, athough, this is largely a moot point because Gaussian ’09 includes Grimme’s dispersion correction for optimisation (for some functionals).

References

Grimme, S. (2006). Semiempirical gga-type density functional constructed with a long-range dispersion correction. Journal of Computational Chemistry, 27(15):1787–1799.

Grimme, S., Ehrlich, S., and Goerigk, L. (2011). Effect of the damping function in dispersion corrected density functional theory. Journal of Computational Chemistry, 32(7):1456–1465.

The bash script:

#!/bin/bash

for file in $*

do
./djag2-3 $file >> output
done

exit 0

The Perl script

 #!/usr/bin/perl

#use Devel::Dprofiler;

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# djag: Batch Dispersion Calculator for Jaguar Output v2.3 #
# ------------------------------------------------ #
# #
# Benjamin M. Ridgway #
# #
# Centre for Computational Chemistry, #
# University of Bristol, 2012. #
# #
# benjamin_ridgway (at) me.com #
# #
# a script to calculate the van der Waals contribution to #
# the energy of a molecule #
# #
# accepts atomic coordinate input from Jaguar output #
# #
# For further details see: #
# - Grimme, S.; J. Comput. Chem.; 25(12), 2004. #
# #
# Modified parameters and function forms from: #
# - Grimme, S.; J. Comput. Chem.; 27(15); 2006. #
# #
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#

#+Version History++++++++++++++++++++++++++++++++++++++++++#
#
# Version 1.0-2.0 (BR)
# Added functions - eval statements to catch division by zero errors
# Added - function to catch X atomic coordinates i.e. dummy atoms
# should skip over these lines and only take true atomic coordinates
#
# Version 2.1: (BR) 19-03-2012 (l235-l244)
# Corrected reading incomplete files which leads to division by zero errors
#
# Version 2.2 (BR) 21-03-2012
# Added some extra checks. Will abort if file has completed but contains no 'final geometry' i.e. it was a sp calculation
# or something other than a geometry optimisation.
#
# Version 2.3 (BR) 23-04-2012
# Will now import Jaguar 'in' files
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#

my $filePath;
my @filePathTemp;
my $filePathTempSize;
my $fileName;
my $fileTypeSize;
my @fileTypeComponents;
my $fileExt;

# PARAMETERS
# From Grimme's 2006 paper:

# General parameters

my $d = 20; # 23 in the original (2004) paper

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# units are J nm^6 mol^-1 - later these numbers are converted to hartrees
my %c6param = (H => "0.14",
He => "0.08",
Li => "1.61",
Be => "1.61",
B => "3.13",
C => "1.75",
N => "1.23",
O => "0.7",
F => "0.75",
Ne => "0.63",
Na => "5.71",
Mg => "5.71",
Al => "10.79",
Si => "9.23",
P => "7.84",
S => "5.57",
Cl => "5.07",
Ar => "4.61",
K => "10.8",
Ca => "10.8",
Sc => "10.8",
Ti => "10.8",
V => "10.8",
Cr => "10.8",
Mn => "10.8",
Fe => "10.8",
Co => "10.8",
Ni => "10.8",
Cu => "10.8",
Zn => "10.8",
Ga => "16.99",
Ge => "17.1",
As => "16.37",
Se => "12.64",
Br => "12.47",
Kr => "12.01",
Rb => "24.67",
Sr => "24.67",
Y => "24.67",
Zr => "24.67",
Nb => "24.67",
Mo => "24.67",
Tc => "24.67",
Ru => "24.67",
Rh => "24.67",
Pd => "24.67",
Ag => "24.67",
Cd => "24.67",
In => "37.32",
Sn => "38.71",
Sb => "38.44",
Te => "31.74",
I => "31.5",
Xe => "29.99",
Au => "29.99",
Pt => "29.99");
#Au => "1409.1",
#Pt => "1409.1");

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#

# the name of the output file, *.out

my $numberOfFiles = @ARGV;

my @files = ;

my $fileCounter;

# for ($fileCounter = 0; $fileCounter "1.001",
He => "1.012",
Li => "0.825",
Be => "1.408",
B => "1.485",
C => "1.452",
N => "1.397",
O => "1.342",
F => "1.287",
Ne => "1.243",
Na => "1.144",
Mg => "1.364",
Al => "1.639",
Si => "1.716",
P => "1.705",
S => "1.683",
Cl => "1.639",
Ar => "1.595",
K => "1.485",
Ca => "1.474",
Sc => "1.562",
Ti => "1.562",
V => "1.562",
Cr => "1.562",
Mn => "1.562",
Fe => "1.562",
Co => "1.562",
Ni => "1.562",
Cu => "1.562",
Zn => "1.562",
Ga => "1.65",
Ge => "1.727",
As => "1.76",
Se => "1.771",
Br => "1.749",
Kr => "1.727",
Rb => "1.628",
Sr => "1.606",
Y => "1.639",
Zr => "1.639",
Nb => "1.639",
Mo => "1.639",
Tc => "1.639",
Ru => "1.639",
Rh => "1.639",
Pd => "1.639",
Ag => "1.639",
Cd => "1.639",
In => "1.672",
Sn => "1.804",
Sb => "1.881",
Te => "1.892",
I => "1892",
Xe => "1.881",
Au => "1.66",
Pt => "1.72");

# counters for finding the final geometry
my $firstMarker;
my $lastMarker;

# some more counters
my $count = 0;
my $i;
my $j; # for creating a new array of size j.

### FROM exJag program

my $jagTemp;
my @newJag;
my @array;

# Remains zero if the job does not contain the expression "completed on"
my $completedSuccessfully = 0;
my $multipleStructuresIsTrue = 0;
my @jagAtomLabel;

THEFILE: foreach (@files) {

$filePath = $_;
@filePathTemp = split(/\//, $filePath);
$filePathTempSize = $filePathTemp;
$fileName = $filePathTemp[$filePathTempSize - 1];
#print "$fileName\n";

@fileTypeComponents = split(/\./, $fileName);
$fileTypeSize = @fileTypeComponents;
$fileExt = $fileTypeComponents[$fileTypeSize - 1];
#print "$fileExt\n";

open( JAGOUT, "<$_" );

@jagout = ; # put each line into an array

close(JAGOUT);

my @tempEnergy;

# COORDINATE EXTRACTION
# if ($jagout[$#jagout] =~ /completed on/){
# print "$files[$#files]: Error reading file, perhaps this calculation ended early or terminated abnormally\n";
# #break;
# continue else;
# my $final_disp = 0;
# my $corE = 0;
# my $array_size = 0;
# exit;
# }

my $readNext = 0;
my @inFileAtomArray;

# ++++++ FOR FILES OF TYPE JAGUAR 'IN'+++++++++

if ($fileExt =~ /in/){

&infile(@jagout); #call the Jaguar in reading routine

}

else {

# ++++ for Jaguar 'OUT' files++++++

foreach (@jagout) {

if ($_ =~ /completed on/){
$completedSuccessfully++;
}

if ($_ =~ /final geometry/){
$multipleStructuresIsTrue++;
}

}

if (($completedSuccessfully == 0) || (($completedSuccessfully != 0) && ($multipleStructuresIsTrue == 0))){

print "$files[$#files]: did not complete successfully.\n";

$completedSuccessfully = 0;

$multipleStructuresIsTrue = 0;

exit;
}

foreach (@jagout) {

if ($_ =~ / SCFE: SCF energy:/) {
@tempEnergy = split( /\s+/, $_ );

}
}

$jaguarEnergy = $tempEnergy[5];

$count = 0;

foreach (@jagout) { # For each of the lines of the file read in.

if ( $_ =~ /^ final geometry/ )
{ # Find the string matching 'final geometry', where

$firstMarker =
$count; # Mark this location using the counter so we can find it later

}

# if ( $_ =~ /^ principal moments of inertia:/ )
# { # Match the very end of the molecule specification...
#
# $lastMarker = $count; # ...and make a note of its location
#
# }

$count = $count + 1; # Increment the counter by one

}

#print "$_, $firstMarker, $lastMarker\n";
# ARRAY PROCESSING
my $atomType;

# some new scalars
my $fileLine = $firstMarker + 3; # the first line of the atomic coordinates

my $fileAtomicCounter = 0; # an array index counter

my $fileAtomic; # holds the file's atomic coordinate lines

my $nextLineBlank = 0;

while ( $nextLineBlank == 0 ) {

#if ($jagout [ $fileLine ] =~ /[A-Z][a-z][0-9]{1,3}\s/ ) { # matches e.g. Pd1 or Pd123

$fileAtomic = $jagout[ $fileLine ]; # put the line in the array

$fileAtomic =~ s/\s\n//g; #remove any whitespace and the trailing newline

$fileAtomic =~ s/^\s{1,2}//g; #remove any leading whitespace

$fileAtomic =~ s/\s+/;/g; # replace the remaining whitespace with a semicolon
# so that we can split the field into an array

@newJag = split( /\;/, $fileAtomic ); # make an array of the data, size 4

$array[$fileAtomicCounter] = [@newJag]; # create a new array, each line containing the previously generated array

$jagAtomLabel[$fileAtomicCounter] = $array[$fileAtomicCounter][0];

$atomType[$fileAtomicCounter] = $array[$fileAtomicCounter][0];

$atomType[$fileAtomicCounter] =~ s/[0-9]+$//g;
#print "$atomType[$fileAtomicCounter]\n";
$fileAtomicCounter++; # increment the array index counter

if ( $jagout [ $fileLine + 1 ] =~ /^\s+\n/ ) { # if the following line is spaces and newline

$nextLineBlank = 1;

} else {

$fileLine++;

}

#}

}

#####

# SCALARS

my $edft=0; # the SCF energy from Gaussian (E(DFT))
my $current_line; # $_
my @iteration_array; # each line matching 'SCF Done' string as array
my @gaulog; # an array that holds the...
# ...contents of the Gaussian log

my $filestem = $ARGV[0]; # the name of the inputfile, the first perl argument
my $log = ".log"; # extension of the logfile
my $file = $filestem.$log; # concatenate the two previous variables

my $outputfile = $filestem."vdw"; # not used in this program (for gaussian output)

# END READ LOGFILE SECTION
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#

}
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#

# SCALARS

my @mDist_ij; # the distance matrix
my @mC_ij; # the Cij coefficient matrix
my @mDamp; # the damping function matrix
my @mRadii_ij; # the van der Waals matrix

my @deriv_temp; # holds the temporary first derivative information
my @matrix_d; # this doesn't seem to be used anywhere
my $T; # another counter

my @forces;

my $i; # a counter for atom i
my $j; # a counter for atom j

# these variables are used to compute
# the distance between atom pairs

my $xyz; # a counter for xyz coordinates
my $diff = 0; # temporary variable holding x(i)-x(j)
my $sqdiff = 0; # holds (x(i)-x(j))^2
my $sumsq = 0; # holds the square root of $sqdiff
my $damp_ij;
my $c_ij;
my $r_ij;
my $pairs; # the number of unique atom pairs

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#

# FILE INPUT SECTION

# SCALARS

# a pair of temporary arrays, these are used to separate
# out the parts of the inputfile

my @temp_array;
my @temp_array2;

my @atom_array; # atom array holds atomic no, x, y, z
my @atom; # atom holds only x, y, z

my @param; # the contents of the parameter file
my @param_array; # the parameter file as an array of arrays
my @atomic_number; # the atomic numbers of the atoms from inputfile
my @c6_array; # the C6 coefficients of the input atoms
my @rij_array; # the van der Waal radii of the input atoms
my @distance_matrix; # matrix contains the dx, dy and dz distances

# some counters

my $t;
my $u;

my $temp = 0 ;

# FILE INPUT FUNCTIONS

my $array_size = @array; # number of lines in the array

my $array2_size = @param;

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#

# GET THE COORDINATES

for ( $u = 0 ; $u <= $array_size-1 ; $u++ ) {

for ( $t = 0 ; $t <= 3 ; $t++ ){

# the current line in the array (one line from input)

# $temp = $array[$u];

# split the line into an array (split at whitespace)

# @temp_array = split(/\s+/, $temp);

# @atom_array contains coords

$atom_array[$u][$t] = $array[$u][$t];

}

}

my $max = @atom_array; # the number of atoms in the molecule

# GET THE PARAMETERS

for ( $u = 0 ; $u <= $array2_size-1 ; $u++ ) {

for ( $t = 0 ; $t <= 2 ; $t++ ){

# as before read each line from the input array

$temp2 = $param[$u];

# and split it into its components
# in this case the atomic number, C6, and vdW radius

@temp2_array = split(/\s+/, $temp2);

# makes an array of arrays of the input parameters
# so that we can look up the parameters for each
# element from the input coordinate file

$param_array[$u][$t] = $temp2_array[$t];

}

}

for ( $u = 0 ; $u <= $max - 1 ; $u++ ) {

for ( $t = 1 ; $t <= 3 ; $t++ ){

# strip off the first element of each array of array
# $atom_array[u][0] is the atomic number
# make a new array of arrays of x, y, z coordinates

$atom[$u][$t-1] = $atom_array[$u][$t];

}
#print "$atom[$u][0], $atom[$u][1], $atom[$u][2]\n";
}

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++#
# GET OTHER PARAMETERS (c6 rij)

# turn the atomic number array into an array of C6 coefficients

my $counter = @array;

for ( $i = 0; $i <= $counter-1; $i++ ) {

$c6_array[$i] = $c6param{$atomType[$i]};

$rij_array[$i] = $vdw{$atomType[$i]};

#print "$atomType[$i] : $c6_array[$i] $rij_array[$i]\n";
}

for ( $i = 0 ; $i <= $max - 2 ; $i++ ) {

for ( $j = $i + 1 ; $j <= $max - 1; $j++ ) {

for ( $xyz = 0 ; $xyz <= 2 ; $xyz++ ) {

# find the difference x(i)-x(j)

$diffxyz =
$atom[$i][$xyz] - $atom[$j][$xyz];

# find the square of $diffsyz

$sqdiff = $diffxyz * $diffxyz;

# add up the x, y and z components
$sumsq =
$sqdiff + $sumsq;
}

# and find the square root - the separation (angstroems)

my $dist = sqrt $sumsq;

#21-03-2012 the following line for killing for division by zero. Probably need to delete this at some point.
#if ($dist == 0) { exit;}
# print "$atomType[$i]-$atomType[$j]: $dist\n";

$mDist_ij[$i][$j] = $dist; # matrix element ij contains
# the squared differences
# for the atom pair should
# be (i*j)/2 elements...
#print "i $i j $j is $mDist_ij[$i][$j]\n";

# reinitialize variables for next atom j

$sqdiff = 0;
$sumsq = 0;
$dist = 0;

}

}
my $mDist_ij_count = @mDist_ij;
#print "mDist_ij contains $mDist_ij_count\n";

# calculate a matrix of C6 coefficients and vdW sums
# iterate over the atom pairs atom(i)atom(i+1) over all i

my $hartree_in_j_per_mol = 2625499.62;

for ( $i = 0 ; $i <= $max - 2 ; $i++ ) {

for ( $j = $i + 1 ; $j <= $max - 1; $j++ ) {

# earlier form of the C6 equation (from 2004)

#$rk = 2 * (( $c[$i] * $c[$j] ) / ( $c[$i] + $c[$j] ));
#$mC_ij[$i][$j] = 1.2 * $rk;

# modified geometric mean form 2006 paper

eval { $c_ij = sqrt( $c6_array[$i]/$hartree_in_j_per_mol * $c6_array[$j]/$hartree_in_j_per_mol );
$mC_ij[$i][$j] = $c_ij;}; warn $@ if $@;

# calculate the sums of the van der Waal radii
# for all the atom pairs

eval {$r_ij = $rij_array[$i] + $rij_array[$j];
$mRadii_ij[$i][$j] = $r_ij;}; warn $@ if $@;

}

}

my $mRadii_ij_count = @mRadii_ij;
my $mC_ij_count = @mC_ij;

for ( $i = 0 ; $i <= $max - 2 ; $i++ ) {

for ( $j = $i + 1 ; $j <= $max - 1; $j++ ) {

eval { $damp_ij =
1 / ( 1 + exp( -$d * ( ( $mDist_ij[$i][$j] /
$mRadii_ij[$i][$j] ) - 1 ) ) ); }; warn $@ if $@;

$mDamp[$i][$j] = $damp_ij;

}

}
my $mDamp_count = @mDamp;
#print "mDamp contains $mDamp_count\n";
# calculate the dispersion correction for all the atoms

# temporary variable for holding one iteration of Edisp

my $edisp = 0;

# iterate over the atom pairs atom(i)atom(i+1) over all atoms

for ( $i = 0 ; $i <= $max - 2 ; $i++ ) {

for ( $j = $i + 1 ; $j <= $max - 1; $j++ ) {

eval { $disp = ( $mC_ij[$i][$j] / ( ($mDist_ij[$i][$j])/10)**6 )
* $mDamp[$i][$j] ; }; warn $@ if $@;

$edisp = $edisp + $disp;

}

}

# calculate the total dispersion correction by multiplying
# by the scaling factor

my $final_disp = $edisp * (-1.05) ;

my $corE = $jaguarEnergy + $final_disp;

my @filenamer = split( /\//, $_ );

print "$filenamer[$#filenamer] $array_size $final_disp $jaguarEnergy $corE\n";
my $final_disp = 0;
my $corE = 0;
my $array_size = 0;
}
# Subroutine INFILE
# returns coordinates from a jaguar 'in' file as array @array with each element
# containing an array of atom label, x, y and z.

sub infile() {

foreach (@_){
if ($_ =~ /&\n/) {

$readNext = 0;
}

elsif ($readNext == 1) {

push(@inFileAtomArray, $_)

}

elsif ($_ =~ /&zmat/) {

$readNext = 1;

}

}

foreach (@inFileAtomArray){
$_ =~ s/\s?\n//g; #remove any whitespace and the trailing newline

$_ =~ s/^\s{1,2}//g; #remove any leading whitespace

$_ =~ s/\s+/;/g; # replace the remaining whitespace with a semicolon
# so that we can split the field into an array

@newJag = split( /\;/, $_ ); # make an array of the data, size 4

$array[$fileAtomicCounter] = [@newJag]; # create a new array, each line containing the previously generated array

$jagAtomLabel[$fileAtomicCounter] = $array[$fileAtomicCounter][0];

$atomType[$fileAtomicCounter] = $array[$fileAtomicCounter][0];

$atomType[$fileAtomicCounter] =~ s/[0-9]+$//g;

$fileAtomicCounter++;

}

return $array;
}
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s