[gmx-users] About the QM/MM MD of Gromacs
Clinton King
clintonking36 at chem.byu.edu
Fri Jul 28 21:13:25 CEST 2017
In case anyone is interested, here is a copy of the script I am currently
using to interface Gromacs 2016.1 with Gaussian 09 without altering the
source code of Gaussian 09. It is based on Gerrit Groenhoef's script,
available at wwwuser.gwdg.de/~ggroenh/qmmm.html#gaussian, with some
alterations so that it better fits my particular system. One note, since I
use custom basis sets with Gaussian 09, this script requires that I include
a file in the job submission with the extension .basis so that the basis
set information is concatenated to the end of the job file that gets
submitted to Gaussian.
The supercomputer at our institution uses slurm for job submission, and I
set the following variables in my slurm job submisstion script:
export OMP_NUM_THREADS=1
export Data_prefix=/fslhome/crk9/gromacs/gromacs-2016.1/bin/gmx_mpi #
location of Gromacs 2016.1 executable
export GMX_QM_GAUSS_EXE=gau
export GMX_QM_GAUSS_DIR=/fslhome/crk9/gromacs/ # the location of gau
export GMX_QM_GAUSSIAN_NCPUS=16
export GMX_QM_GAUSSIAN_MEMORY=2 # part of a workaround, since I couldn't
figure out how to tell gromacs that I wanted to give gaussian 16 GB of
memory to use
export QM_METHOD=M06
export QM_MEM=16gb
export QM_PROC=16
*start of script:*
#!/bin/bash
#
# Please read this carefully, you may need to change it on a few locations
#
# wrapper to call gaussian, convert in- and output with some (tedious)
# awk scripting into the 'old' formats and perform the call to the
# standard g09/g03 installed on the system. For reasons I don't
# understand the system() routine does not like long filenames. I
# therefore call this script gau
#
# To make gromacs aware of this script, use GAUSS_EXE:
#
# export GMX_QM_GAUSS_EXE=
#
# Also make sure that gromacs can find it by using the full path
# location to this file:
#
# export GMX_QM_GAUSS_DIR=
#
#the next three lines are to include variables to fix the qm method and the
memory for gaussian, also a log file is generated to help track what's
going on in the qm #computations
echo "qm method is $QM_METHOD" >> QM.stats.log
echo "amount of memory for gaussian is $QM_MEM" >> QM.stats.log
echo "number of processors is $QM_PROC" >> QM.stats.log
# let gromacs write a "normal" gaussian input file, using the keyword
# determined by the User's input. Alternatively, it is possible now to
# provide different keywords by placing a header.com in the
# directory. This will take precedence over the gromacs keywords.
echo
'*********************************************************************' >>
QM.stats.log
echo 'new run' >> QM.stats.log
if [ -e QM.step ]
then
step=$(head -n 1 QM.step)
#echo "$step =>" >> QM.stats.log
[[ "$step" =~ ([0-9]+)$ ]] && step="$((BASH_REMATCH[1] + 1))";
#echo "$step" >> QM.stats.log
echo "Current step = $step" >> QM.stats.log
else
step=0
fi
echo $step > QM.step
# create some temporary files
if [ -e temp.com ]
then
rm -fr temp.com
echo 'removed temp file' >> QM.stats.log
else
echo 'temp file does not exist' >> QM.stats.log
fi
# The user has provide an own header.com with keywords, including the
title, but no whiteline after that!
if [ -e "$JOBNAME".com ]
then
cat "$JOBNAME".com >> temp.com
echo 'jobname = '>> "$JOBNAME" >> QM.stats.log
echo 'jobname.com exists, copied it to temp.com' >> QM.stats.log
else
# Use the gromacs keywords based on input in the mdp file
echo '--------input.com-----------'>>QM.stats.log
cat input.com >> QM.stats.log
echo '------------'>>QM.stats.log
echo "%nprocshared=$QM_PROC" >> temp.com
awk -v w=0 '{if ( ( w == 0 ) && ( $1 != "%subst" )) print}{ if ($1 == "#P"
|| $1 == "#T") w=1}' input.com >> temp.com
echo "Nosymm units=bohr FORCE Punch=(Derivatives) iop(3/33=1)
Prop=(Field,Read) pseudo geom=nocrowd" >> temp.com
echo " " >> temp.com
echo "empty title" >> temp.com
echo "jobname = $JOBNAME" >> QM.stats.log
echo 'jobname.com exists, copied it to temp.com' >> QM.stats.log
fi
# Use the knowledge on the input structure to extract the coordinates, etc.
awk -v w=0 '{if (w == 1) print}{if ( $1 == "input-file") w=1}' input.com >>
temp.com
echo 'coordinates extracted' >> QM.stats.log
# The following section will add basis set information necessary when using
the 'gen' keyword
if [ -e "$JOBNAME".basis ]
then
cat "$JOBNAME".basis >> temp.com
echo 'basis exists and copied to temp.com' >> QM.stats.log
else
echo 'basis does not exist' >> QM.stats.log
fi
# IN order to get the gradients on the MM point charges, we request
# gaussian to compute the electrostatic gradients on these centers. We
# thus need to provide the coordinates of the nearby MM charges once more.
We
# again use our knowledge of the input structure and simply wait until
# we have read in four whitelines, which means we have arrived at the
# point charges field.
# Note that gaussian ALWAYS requires these positions to be in angstrom,
# even if units=bohr!
awk -v nlines=0 '{if (nlines == 4 && $1) printf "%10.8f %10.8f %10.8f\n",
$1*0.529177249,$2*0.529177249,$3*0.529177249}{if ($1) ;else nlines++ }'
input.com >> temp.com
echo 'point charges field extracted' >> QM.stats.log
#adds a blank line to the end of gaussian input file
echo "" >> temp.com
# We need another temporary file to remember the charges, which we
# will multiply later when the gaussian call is done with the
# potential to get the gradients. (These are the charges of the QM nuclei
and nearby MM atoms)
if [ -e MMcharges.txt ]
then
rm -rf MMcharges.txt
echo 'MM charges did exist' >> QM.stats.log
else
echo 'MM charges did not exist' >> QM.stats.log
fi
touch MMcharges.txt
echo "touched MMcharges" >> QM.stats.log
awk '{if ($1 !~ /\./ ) {if ($1 ~ /^[0-9]/) {if ($4) print $1;};} else
print $4}' input.com >> MMcharges.txt
echo "copied input to charges" >> QM.stats.log
sed -i -e "s/DFT/$QM_METHOD/g" temp.com
echo "replaced DFT with $QM_METHOD"
#always set QMbasis in mdp file to STO-3G
sed -i -e 's/STO-3G/gen/g' temp.com
echo "replaced STO-3G with gen"
#This is a workaround for communicating memory to gaussian (always set
GMX_QM_GAUSSIAN_MEMORY=2)
sed -i -e "s/%mem=2/%mem=$QM_MEM/g" temp.com
echo "set mem to $QM_MEM for Gaussian"
# Now the input file is in principle ready and we can make a call to
# gaussian. If additional input fields are required after the
# pointcharges, such as the state averaging coefficients, they can be
# entered here via another cat command.
echo "READY TO SUBMIT JOB TO GAUSSIAN" >> QM.stats.log
cat temp.com >> QM.stats.log
/fslapps/chem/bin/rung09 temp.com
export return=$status
# build in a check if gaussian quit cleanly (look for the phrase 'Normal
termination' on last line,
#or 'Error termination' on third to last line.)
# after finishing, we can read in the input
echo "FINISHED GAUSSIAN RUN" >> QM.stats.log
if [ -e temp.7 ]
then
rm -rf temp.7
echo 'temp.7 existed, and removed' >> QM.stats.log
fi
touch temp.7
echo 'temp.7 touched'>>QM.stats.log
#self-energy of MM charges need to be subtracted from total energy
grep "Self energy" temp.log |awk '{print $7}' > pointchargeself.txt
echo 'grabbing self energy from temp.log' >> QM.stats.log
# look for the final energy
grep "SCF Done" temp.log | awk '{print $5}' > SCF.txt
echo 'grabbing final energy from temp.log' >> QM.stats.log
paste SCF.txt pointchargeself.txt|awk '{print $1-$2}' >> temp.7
echo 'paste energies to temp.7' >> QM.stats.log
# gradients on the QM atoms
cat fort.7 >> temp.7
echo 'cat fort.7 (contains gradients on QM atoms) to temp.7' >> QM.stats.log
#get the gradients of the QM atoms and nearby MM atoms
awk -v start=0 '{if ($5 == "Field") start=NR+3}{if(start==NR) start=1}{if
(start==1 && NF==1) start=0}{if (start==1 && NF>=5) print}' temp.log >
gradients.txt
#Gaussian outputs the electric field in atomic units, and gromacs needs the
electric field
#to be in kJ mol-1 nm-1 e-1. The numbers from Gaussian need to be
multiplied by 2.0155295489e-5, but this is taken care of in the gromacs code
paste MMcharges.txt gradients.txt|awk '{if ($3 == "Atom") print
-1*$1*$5,-1*$1*$6,-1*$1*$7; else print -1*$1*$4,-1*$1*$5,-1*$1*$6}' |
column -t >> temp.7
#prints the fort.7 before we replace it later
#cat fort.7 > $JOBNAME-preFORT-$step.7
rm -rf fort.7
echo 'gets the energies of MM atoms (and other things)' >> QM.stats.log
# Guassian uses D rather than e in scientific representation, so we change
that.
sed 's/D/E/g' temp.7 > fort.7
cat temp.com
cat temp.log
cat fort.7
cat temp.7
#prints the submission file and result for user reference.
#cat SCF.txt > $JOBNAME-SCF-$step.txt
#cat pointchargeself.txt > $JOBNAME-pointchargeself-$step.txt
#cat gradients.txt > $JOBNAME-gradients-$step.txt
#cat MMcharges.txt > $JOBNAME-MMcharges-$step.txt
#cat $JOBNAME.basis > $JOBNAME-basis-$step.txt
#cat input.com > $JOBNAME-GRO-$step.com
#cat temp.com > $JOBNAME-QM-$step.com
#cat temp.log > $JOBNAME-QM-$step.log
#cat fort.7 > $JOBNAME-FORT-$step.7
exit $return
*end of script*
--
Clinton King
Graduate Student
Chemistry Department
Brigham Young University
On Fri, Sep 9, 2016 at 10:03 AM, Clinton King <clintonking36 at chem.byu.edu>
wrote:
> An FYI: I have found some bugs in the awk scripting on (about) lines 59
> and 71 of the gau script. The following two lines of awk code should fix
> those bugs, but as soon as I am sure I have it working, I'll post the
> version of the gau script that I am using to this forum in case anyone else
> wants to use it.
>
> awk -v nlines=0 '{if (nlines > 2 && $1) printf "%10.8f %10.8f %10.8f\n",
> $2*0.529177249,$3*0.529177249,$4*0.529177249}{if ($1) ;else nlines++ }'
> input.com >> temp.com
>
> awk -v nlines=0 '{if (nlines > 2 && $1) print $1}{if ($1) ;else nlines++
> }' input.com >> MMcharges.txt
>
> --
> Clinton King
> Graduate Student
> Chemistry Department
> Brigham Young University
>
>
>>
>> Today's Topics:
>>
>> 1. Re: About the QM/MM MD of Gromacs (???) (Timofey Tyugashev)
>>
>>
>> ----------------------------------------------------------------------
>>
>> Message: 1
>> Date: Fri, 9 Sep 2016 15:50:49 +0700
>> From: Timofey Tyugashev <tyugashev at niboch.nsc.ru>
>> To: gromacs.org_gmx-users at maillist.sys.kth.se
>> Subject: Re: [gmx-users] About the QM/MM MD of Gromacs (???)
>> Message-ID: <73e2df17-6226-35bd-688b-c1413bfce167 at niboch.nsc.ru>
>> Content-Type: text/plain; charset=UTF-8; format=flowed
>>
>>
>> So the only QM software that GROMACS can run for QMMM is Gaussian (09?),
>> in single-thread CPU-only mode and with cut-off setting that's listed as
>> depreciated, right?
>> 08.09.2016 21:11, gromacs.org_gmx-users-request at maillist.sys.kth.se
>> ?????:
>> > Message: 3
>> > Date: Thu, 8 Sep 2016 13:05:38 +0000
>> > From: "Groenhof, Gerrit"<ggroenh at gwdg.de>
>> > To:"gromacs.org_gmx-users at maillist.sys.kth.se"
>> > <gromacs.org_gmx-users at maillist.sys.kth.se>
>> > Subject: Re: [gmx-users] About the QM/MM MD of Gromacs (???)
>> > Message-ID:
>> > <858A7947BC0FE04DA05E1786A6D51D4533080324 at um-excdag-a05.um.
>> gwdg.de>
>> > Content-Type: text/plain; charset="us-ascii"
>> >
>> > Hi,
>> >
>> > Orca with gromacs >= 5 has not been tested and there are some
>> indications that it is broken (redmine issue 1934).
>> >
>> > QM/MM with Gaussian is working, as long as you use a group based
>> cut-offs and run mdrun on a single thread (-nt 1).
>> >
>> > I therefore recommend to use gaussian (with the gau script athttp://
>> wwwuser.gwdg.de/~ggroenh/qmmm.html#gaussian no source code of Gaussian
>> is required), or stick to an older version of gromacs that supports ORCA.
>> >
>> > Best,
>> >
>> > Gerrit
>> >
>> >
>> >
>> >
>> > Dear all:
>> > I want to do the QM/MM MD calculation by gromacs. ORCA is my QM
>> program for QM/MM calculation. But how can I configure this option by cmake?
>> > For older version of gromacs, we just need to run the command:
>> > ./configure --with-qmmm-orca --without-qmmm-gaussian
>> > But for newer version of gromacs, we can only use the cmake to
>> configure this option. Then how can I do? What is the command?
>> > Thank you!
>> > Best wishes,
>> > Junbo
>>
>>
>>
>>
>>
>
More information about the gromacs.org_gmx-users
mailing list