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 <clintonkin...@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 <tyugas...@niboch.nsc.ru> >> To: gromacs.org_gmx-users@maillist.sys.kth.se >> Subject: Re: [gmx-users] About the QM/MM MD of Gromacs (???) >> Message-ID: <73e2df17-6226-35bd-688b-c1413bfce...@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-requ...@maillist.sys.kth.se >> ?????: >> > Message: 3 >> > Date: Thu, 8 Sep 2016 13:05:38 +0000 >> > From: "Groenhof, Gerrit"<ggro...@gwdg.de> >> > To:"gromacs.org_gmx-users@maillist.sys.kth.se" >> > <gromacs.org_gmx-users@maillist.sys.kth.se> >> > Subject: Re: [gmx-users] About the QM/MM MD of Gromacs (???) >> > Message-ID: >> > <858a7947bc0fe04da05e1786a6d51d4533080...@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 >> >> >> >> >> > -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.