This is flucq, produced by makeinfo version 4.0 from flucq.texi.

Combined QM/MM Fluctuating Charge Potential for CHARMM Ben Webb, ben@bellatrix.pcl.ox.ac.uk, and Paul Lyne The fluctuating charge potential (FlucQ or FQ) is based on the method developed by Rick, Stuart and Berne (Rick et. al., J. Chem. Phys. 101 (7) 1994 p6141) for molecular dynamics, and extended for hybrid QM/MM simulations (Bryce et. al., Chem. Phys. Lett. 279 1997, p367). It is designed primarily for computationally efficient (approx. 10% overhead) modelling of solvent polarisation in hybrid QM/MM systems, and as such is implemented for QUANTUM, CADPAC and GAMESS codes, although the current implementation is easily extensible to any atom type and bond. * Menu: * Syntax:: Syntax of the FLUCQ command * Activation:: Starting FlucQ from a CHARMM input file * Charge solution:: Solving for exact charges * Reference energy:: Setting the ``zero'' for FlucQ polarisation * Caveats:: Changes to be aware of; known limitations * Using FlucQ with QM:: Necessary changes for use with CADPAC or GAMESS * Examples:: Simple uses of the FLUCQ command * Implementation:: Mathematical and computational details

[SYNTAX FLUCq] FLUCq { ON init-spec (atom selection) } { OFF } { PRINt } { EXACt exac-spec } { REFErence { GAS exac-spec } } { { SOLVent exac-spec } } { { CURRent } } { { ENERgy real } } init-spec::= [GROUp] [NOFIxed] nose-spec nose-spec::= [TEMPerature real] [TMASs real] [TOLEranace real] [ITERations integer] exac-spec::= [TIMEstep real] [ZETA real] [TQDEsired real] [PRINt]

FlucQ code is enabled within CHARMM by means of the FLUCQ ON command. Future energy calculations will then include an extra energy term - FQPO, the FlucQ polarisation energy, while dynamics simulations involve a new energy property - FQKI, the FlucQ charge kinetic energy. Once FlucQ is active, the selected atoms are treated as extra degrees of freedom, free to fluctuate under the charge forces in the system, and, by assigning each atom type a fictional charge "mass", these charges can be accelerated in a conventional dynamics simulation, in a completely analogous way to the Cartesian degrees of freedom. If atoms are selected by the FLUCQ command which cannot be modelled (i.e. they are QM atoms, or have no FlucQ parameters defined for them) they will be automatically removed from the selection. The FlucQ polarisation energy, FQPO, is an intramolecular interaction; in full electronegativity equalisation, every atom interacts through space, by means of a modified Coulomb-type interaction, with every other atom in the molecule. In this implementation, the only interactions calculated are those along defined CHARMM bonds (even those with zero force constants). [TEMPerature <real>] specifies the temperature at which charge degrees of freedom will be maintained during dynamics by Nose-Hoover thermostatting (default 0). The fictional mass associated with the thermostat can be set with the TMASs option (default 0, no thermostatting), the tolerance of the Nose-Hoover iterations with TOLErance (default 1.0d-7), and the maximum number of iterations with ITERations (default 100). [GROUp] conserves charge within groups, rather than the default behaviour of conserving charge within residues; this prohibits charge transfer between groups. Note that the FlucQ model makes no restriction on the degree of charge transfer within each residue or group, or the distance over which this transfer can occur. [NOFIxed] instructs FlucQ that some or all of the bond lengths between FlucQ-selected atoms are free to change during a simulation. This forces the FlucQ code to recalculate the intramolecular interaction at each step; since this is a costly calculation, the default is to use interactions parameterised for equilibrium bond lengths, with which it is strongly recommended to combine constraint methods such as SHAKE BONH PARA. The FLUCq PRINt command simply prints the current values of all charges and charge forces (from the last energy calculation). A similar effect can also be achieved with the standard SCALAR command (see scalar.doc for information on other FlucQ parameters available with the SCALAR command). The FLUCQ OFF command disables the FlucQ code. Further energy calculations will not include FlucQ terms. Note, however, that if the charges have been modified by FlucQ, they will remain at their altered values. The initialization process dimensions FlucQ with the current state of the system. The QM region, if any, is detected, and the FlucQ atom selection will then interact with the QM region. Thus, the FLUCQ command should be placed after any QUANTUM, CADPAC, or GAMESS command, and if the total number of atoms in the system is modified, FlucQ should be disabled prior to this change and reinitialized afterwards. To skip FlucQ energy calculations entirely, use the SKIP FQPOL FQKIN command. The QM/MM FlucQ interaction is calculated in line with the standard QM/MM electrostatic interaction, and as such is suppressed with the SKIP QMEL command. Finally, the intermolecular contribution to FlucQ is calculated in line with the standard electrostatic interaction, and so is disabled with the SKIP ELEC command. No FlucQ interaction energies are calculated between atoms constrained with the CONS FIX command, as electrostatic energies are not calculated for these atoms. FlucQ parameters are specified in the parameter file, with the FLUCQ keyword. The section should look like the following:- FLUCQ atom chi zeta prin mass Here, chi is an electronegativity measure (in Kcal/mol/e), zeta a Slater orbital exponent (in 1/Angstrom), prin the Slater orbital principal quantum number, and mass the charge mass (in (ps/e)**2 Kcal/mol) from the FlucQ model. For example, Rick's original parameters for TIP4P hydrogen and M-site would be written as:- FLUCQ HP 10.00 0.90 1 6.0d-5 MP 78.49 1.63 2 6.0d-5

The FlucQ model relies on keeping charge kinetic energy at a temperature close to zero Kelvin, to maintain Born-Oppenheimer separation between it and the other degrees of freedom. Thus, it is best to acquire a minimum energy charge configuration for your system before any dynamics simulation. Two methods are available for such "charge solution". The first is to use a standard CHARMM minimisation; FlucQ charges will be minimised concurrently with the Cartesian coordinates. The second method is to apply dissipative Langevin dynamics to the charges only, to achieve minimum energy charges for fixed atomic coordinates; this is performed by means of the FLUCq EXACt command. The code prints a running count of the number of iterations required to quench the kinetic energy. [TIMEstep real] sets the timestep to be used in Langevin dynamics, by default 0.001ps. [ZETA real] sets the frictional coefficient, by default 1600. [TQDEsired real] sets the desired final temperature, by default 1.0d-6 K. [PRINt] if set, prints the final charges.

By default, the charge polarisation energy FQPOL reported by FlucQ is given relative to all atomic charges being zero. More generally, it is useful to define this term relative to an arbitrary zero. This reference energy can be set with the FLUCQ REFErence command. FLUCQ REFE GAS disables all intermolecular interactions, solves for exact charges, and then uses the resultant energy as the reference. This essentially defines the polarisation energy relative to the energy that the system would have in the gas phase, with all residues or groups infinitely separated. FLUCQ REFE SOLVENT merely disables the QM/MM interaction, and then sets the reference energy similarly. This shows polarisation as a function purely of the QM system. FLUCQ REFE CURRent defines the current polarisation energy (from the last energy calculation) to be zero - i.e. the reference energy is increased by the current energy. FLUCQ REFE ENERgy real sets the reference energy to a user-specified value. Bear in mind that REFE GAS exac-spec is essentially identical to the series of CHARMM commands:- FLUCQ REFER ENER 0 SKIP ALL EXCL FQPOL BOND ANGL UREY DIHE IMPR FLUCQ EXACT exac-spec FLUCQ REFER ENER ?FQPO SKIP EXCL ALL (The only difference is that any SKIP command in force before REFE GAS will remain in force afterwards, whereas the above example will re-include calculation of all energy terms at completion. Also, by changing the second line in the above example to SKIP QMEL QMVDW, the action of the REFE SOLVENT command can be reproduced.)

The fluctuating charge code alters the atomic charges during dynamics runs. Thus, the charges cannot be treated as constant and restart and trajectory files must include atomic charges. Files read or written during FlucQ-enabled dynamics runs will be assumed to contain charge information, and so will be a) somewhat larger and b) incompatible with non-FQ files. (If FlucQ is compiled in but not activated with FLUCQ ON, the restart and trajectory file formats are unchanged from standard CHARMM.) The FlucQ model is implemented primarily for the study of QM/MM systems, with a fluctuating charge SHAKE-constrained MM solvent. Hence, intramolecular interactions are restricted to those between FlucQ atoms along fixed-length bonds. This complicates the application of the model to large systems, as for full electronegativity equalisation, every atom must interact with every other atom in the group. FlucQ is not implemented for all nonbond routines, in particular the CFF, MMFF, CRAYVEC and PARVECT codes. FlucQ with Ewald is also not currently supported.

In order for the QM/MM calculation to be properly calculated, FlucQ requires data to be passed back to it from the QM codes (in particular the density matrix and one-electron integrals). Changes have been made to the QUANTUM interface for this to be carried out correctly; however, the GAMESS and CADPAC codes, not being distributed with CHARMM, may require modification. These modifications will not affect the functioning of standard QM/MM calculations, when FlucQ is disabled. GAMESS requires the addition of one line of code to grd1.src, in the STVDER subroutine. Specifically, the code near line 1270 should now read:- C ----- NUCLEAR REPULSION FORCE ----- C ----- DENSITY FORCE ----- C ----- HELLMANN-FEYNMAN FORCE ----- C ----- INTEGRAL FORCE (AO DERIVATIVE CONTRIBUTION) ----- C CALL VNNDER(X(LDRG),NAT,X(LEF3),NFRPTS) CALL SDER(X(LEPS),L1,L2) IF (NTMO.GT.0) CALL QMADD(X(LDA),X(LEPS),L2,79) CALL HELFEY(X(LDA),L2) CALL TVDER(X(LDA),L2) CALL DENDD1(X(LDA),X(LDB),L2) C Call to FlucQ CALL FQQMMM(X(LDA)) C IF(SOME .AND. * (IECP.GT.0 .OR. IZRF.GT.0 .OR. IPOT.GT.0 .OR. EFLDL)) THEN CALL TSECND(T1) TG = T1-T0 WRITE(IW,9020) TG T0 = T1 END IF For CADPAC, main.F (subroutine ENRGY) should be changed around line 1250 as below:- C C READ ENRGY FROM DUMPFILE C CALL READ3(EN,LDA(ISEX(13)),ISEC13,IFILD) #ifdef FLUCQ C Call to FlucQ CALL FQQMMM(Q) #endif E=ETOT RETURN END C---------------------------------------------------------- SUBROUTINE GRADY(Q,IQ,MAXQ,ETOT,EG) Subroutine ENUC0 in utilxx.F should be completely replaced with the following:- FUNCTION ENUC0(NAT,ZAN,C) IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON/PAR/node,nnodes,maswrk LOGICAL maswrk DIMENSION ZAN(NAT),C(3,NAT) PARAMETER (MAXLAT=25120) C C LATTICE PARAMETERS C COMMON/LATTIC/CLAT(3,MAXLAT),ZLAT(MAXLAT),NLAT,NL2, 1 BREPL(MAXLAT),AREPL(MAXLAT) C C C LOGICAL LFIELD,FIXED,LEX,LDAM12,LDAM13,LDAM23,LDIIS COMMON/SCFBLK/EN,ETOT,EHF,SH1(2),SH2(2),GAP1(2),GAP2(2), 1 D12,D13,D23,CANA,CANB,CANC,FX,FY,FZ, 2 LFIELD,FIXED,LEX,LDAM12,LDAM13,LDAM23,LDIIS, 3 NCYC,ISCHM,LOCK,MAXIT,NCONV,NPUNCH,LOKCYC, 4 FGXX,FGXY,FGXZ,FGYY,FGYZ,FGZZ #ifdef FLUCQ INTEGER CHMIN #endif C DATA ZERO /0.0D0/ C ENUC0=ZERO IF(NAT.EQ.1) GOTO 20 C C THIS IS THE STANDARD NUCLEAR REPULSION TERM C DO 10 I=2,NAT NI=I-1 DO 10 J=1,NI AX=C(1,I)-C(1,J) AY=C(2,I)-C(2,J) AZ=C(3,I)-C(3,J) RR=AX*AX+AY*AY+AZ*AZ 10 ENUC0=ENUC0+ZAN(I)*ZAN(J)/DSQRT(RR) C C 20 IF(.NOT.LFIELD)GOTO 40 C C THIS PART ADDS INTERACTION OF NUCLEAR CHARGES WITH C EXTERNAL FIELD TO KEEP TOTAL ENERGY CORRECT DURING C FINITE FIELD CALCULATION. C DO 30 I=1,NAT ENUC0=ENUC0-ZAN(I)*(FX*C(1,I)+FY*C(2,I)+FZ*C(3,I)) 30 CONTINUE SXX=0D0 SYY=0D0 SZZ=0D0 SXY=0D0 SXZ=0D0 SYZ=0D0 DO 71 I=1,NAT SXX=SXX+ZAN(I)*C(1,I)*C(1,I) SYY=SYY+ZAN(I)*C(2,I)*C(2,I) SZZ=SZZ+ZAN(I)*C(3,I)*C(3,I) SXY=SXY+ZAN(I)*C(1,I)*C(2,I) SXZ=SXZ+ZAN(I)*C(1,I)*C(3,I) SYZ=SYZ+ZAN(I)*C(2,I)*C(3,I) 71 CONTINUE AV=0.5D0*(SXX+SYY+SZZ) SXX=1.5D0*SXX-AV SYY=1.5D0*SYY-AV SZZ=1.5D0*SZZ-AV SXY=1.5D0*SXY SXY=1.5D0*SXY SXZ=1.5D0*SXZ SYZ=1.5D0*SYZ ENUC0=ENUC0-FGXX*SXX-FGYY*SYY-FGZZ*SZZ 1 -FGXY*SXY-FGXZ*SXZ-FGYZ*SYZ 40 IF(NLAT.EQ.0)RETURN C C If using LATTICE option this adds interaction C of lattice points with each of nuclei in molecule, C but not interactions between lattice points, to keep C total energy correct. C EACH LATTICE POINT HAS A CHARGE (ZLAT) C AND ALSO REPELS EACH REAL NUCLEUS C WITH FORM BREPL * EXP ( -AREPL * R) C DO 60 J=1,NAT #ifdef FLUCQ CHMIN=0 #endif DO 50 I=1,NLAT AX=CLAT(1,I)-C(1,J) AY=CLAT(2,I)-C(2,J) AZ=CLAT(3,I)-C(3,J) RR=DSQRT(AX*AX+AY*AY+AZ*AZ) ENUC0=ENUC0+ZAN(J)*ZLAT(I)/RR 1 + BREPL(I)*DEXP(-AREPL(I)*RR) #ifdef FLUCQ C FlucQ core term CALL FQQCOR(CHMIN,ZAN(J)*ZLAT(I)/RR & +BREPL(I)*DEXP(-AREPL(I)*RR)) C #endif 50 CONTINUE 60 CONTINUE RETURN END

The following example initialises the FlucQ code for a system of SPC waters, before calculating the gas phase energy, and then calculating the self-polarisation of the solvent. Finally, the total energy, including the self-polarisation relative to the gas phase, is printed, and the charge forces from this energy calculation are displayed. FLUCQ ON CHMASS 6.9d-5 TEMP 5.0 SELE RESN SPC END FLUCQ REFER GAS FLUCQ EXACT ENERGY SCALAR FQCFOR SHOW SELE ALL END

The standard CHARMM nonbond routines and QM codes have been modified so as to sum the interaction electrostatic interaction energy between charge "I" and all other nonbond pairs or QM atoms into index "I" of the fluctuating charge array FQCFOR. The FlucQ model actually requires the term dE/dQ, so these totals are divided by charge by the FlucQ energy routine (as all such interactions are linear in charge). Note that this gives erroneous results for FlucQ sites with exactly zero charge; however, the CHARMM nonbond routines calculate no interactions for such systems anyway. Finally, the intramolecular terms, as contributions to dE/dQ, are summed into the FQCFOR array, and charge forces are calculated from these electronegativities by mass-weighted averaging over residues or groups. These forces are then used by the standard minimisers, or by a standard Verlet integrator during dynamics. For further information, see the following:- MM system; (Rick et. al., J. Chem. Phys. 101 (7) 1994 p6141) QM/MM interaction; (Bryce et. al., Chem. Phys. Lett. 279 1997, p367) Tag Table: Node: Top103 Node: Syntax1371 Node: Activation2114 Node: Charge solution6634 Node: Reference energy7793 Node: Caveats9466 Node: Using FlucQ with QM10631 Node: Examples15843 Node: Implementation16424 End Tag Table

Information and HTML Formatting Courtesy of:

NIH/DCRT/Laboratory for Structural Biology

FDA/CBER/OVRR Biophysics Laboratory

Modified, updated and generalized by C.L. Brooks, III

The Scripps Research Institute