ComQum
Comqum is a method to do combined QM/MM calculations.
In principle, it is a general procedure for any QM and MM programs. In practice, it is currently implemented for Amber and Turbomole.

The proper references of ComQum are:

The current version is ComQum-01.
The interface programs are located in:

The code is in /home/bio/Comqum.

ComQum is a free interface program that you can get from Ryde by request.
However, note that to run ComQum, you also need the softwares Turbomole and Amber which are commercial.

The only available documentation is this file, the ComQum page and a detailed specification of the ComQum file formats: comqum01.pdf.



How to start a ComQum calculation (new Amber >6 version)
  Form to fill in
 
  1. Start with an equilibrated crystal structure (with hydrogen atoms) stored in pdb format with charges (follow the instructions on this page).

    So far, we have always run ComQum with sphaerical systems without periodic boundary conditions.
    In principle, it should probably be possible to run with periodic boundary conditions, but this is full of technical problems.
    You must start with pdb structure with the QM system in the centre of the periodic box and all molecules wrapped into the central box (otherwise the point charges will be strange).
    Check also that tleap does not move the atoms (run step 8 on the input pdb file).
    After this, the set-up is normally without problems, but you need to change sander.in2 and 3 to employ periodic boundary conditions and reduce the cutoff to ~12 :
    However, sander often gives infinite energies, owing to problems with the wrapping.
    It is also somewhat unsatisfactorily that the MM calculations are periodic with a finite cutoff, but not the QM calculations.

    If you already have a periodic octahedral simulation, you can truncate it to a spherical system:
    cpptraj prmtop ptraj.in
    with the input file (ptraj.in):
    trajin mdrest3
    trajout mdrest3-wrap restart
    image origin center familiar com :98
    go

    You can get the central residue (98 in this example) by running changepdb, command cap on the pdbout file. First remove all water molecules with command rw and a radius of 0 .

    changeparm <<EOF
    prmtop
    p
    pdb3
    m
    mdrest3-wrap
    q
    EOF

    and then cut it to a spherical system with the cut command of changepdb (use cap first to decide a proper radius, you should include the whole protein and get a spherical system of water molecules, typically 10-15 smaller than the cap radius)
    changepdb
    pdb3
    m
    p
    cu

    0 0 0
    45
    n
    pdb4-cut
    n
  2.  
  3. Change the name of metal-bound water molecules to WAM and move them before the other water molecules. Normally already done

  4. When selecting the quantum system, note that if there is a hydrogen bond from the back-bone NH or CO groups of a junction residue to the quantum system, this group should be included in the quantum system (because they otherwise get strange charges).

  5. Normally not needed: Ensure that you have the desired junctions in the file juncfactor
    (cp /home/bio/Data/Comqum/junctfactor .)

    For each junction, it should contain the following entries:
    Residue name, junction atom, junction neighbour, desired Amber atom type of junction atom;
    Next row: jtype (indicating method and basis set) and ideal bond lenght with this method
    The format is free, but fields must be delimited by spaces.

    $version

      2006

    ASP CB   CA   HC  # Asp OOC-CH3 junction, ff94
    7   1.1065d0      # Acetate, BPRi/6-31G*, Cs, C-H in plane


    Note that the entries depend on the Amber force field to use (the desired atom type) and the QM method to use (the ideal bond length).

    Obsolete: use instead the file /home/bio/Data/Comqum/junctfactor_old
    ASP CB   CA   ACA  1.5260d0 # CT-CT in parm99.dat
     7   1.1065d0               # Acetate, BPRi/6-31G*, Cs, C-H in plane


    Note that this file is specific for a certain force field (e.g. Amber ff99), although this is not explicit.


  6. Run pdbtocomqum
    1. Select a logfile
    2. Enter the name of the pdb file
    3. Do not search for short contatcts
    4. Enter a new title (can be done in a file, see next point).
    5. Select the quantum system, including the junction atoms that will be changed to H atoms (you first give the number of the amino acid, then the number of the atoms in this amino acid; one <CR> let you select a new amino acid; another <CR> ends the selection).
    6. This can optionally (and best) be done by entering the name of a file with the atom numbers (sample file). 
    7. Select the radius for system 2 (MM system that is optimised). Typically ~6 .
    8. Note that it is absolutely essential that system 2 contains the same residues in all QM/MM structures if energies should be comparable. Therefore, you may add or remove residues from the originally suggested list.
    9. Select the radius for system 3 (MM system that is fixed). Typically the rest of the protein.
    10. Select amber (default)
    11. Select the basis set (i.e. the type of junction, typically 12 (TPSS/def2-SV(P)).
    12. Use the new format (default).
    13. Enter the junction atoms (those that will become hydrogen atoms in the quantum system). This can also be done in the file mentioned in point 5.
    14. Check that the junctions are not unknown, i.e. that the junction type is not -1 at the end (if so, look at point 4 above how to insert a new fragment into the libraries).
    15. Do not remove the charge of neighbours (answer n).
    16. Do not redistribute any deleted charges (answer n).

      If you nevertheless decided to remove charges of the neighbours:
    17. Check that the net charge outside system 1 is correct (for normal side-chains, it should be 0, but for backbone fragments and strange groups, it may be <>0).
    18. Check that the change in the charges is not too large. If so, remove more neighbours. For example, you may want to remove charges of hydrogens if you have already removed the charge of their parent heavy atom.
     
  7. Run comqumtoturbo and comqumtoamber
  8. This will give you the following files:
    comqum.dat    comqum.qcout  coord          logfil   pointcharges  sander.in2
    comqum.mmout  control       leap.in        pdb.in3  sander.in1    sander.in3

  9. Check the file coord - it often gives problems with metals (not zn, fe, ni).

  10. Edit the file leap.in to insert unusual residues (i.e. insert extra addPrep and addAmberParm statements) and check the force field. Sample leap.in file.
    Normally, it is best to use the same leap.in file as used for the equilibration run, but you must remove the solvateCap command, change the input pdb file to pdb.in3, and change the output files to prmtop3 and prmcrd3:
    x=loadpdb pdb.in3
    saveamberparm x prmtop3 prmcrd3

    Then, run tleap:

    tleap -f leap.in

    Check that you do not get any errors or warnings, especially after the commands
    x=loadpdb pdb.in3
    saveamberparm x prmtop3 prmcrd3

    Many errors come from missing or extra TER rows in the pdb file - leap will connect all consecutive residues unless they are separated by a TER row.
  11.  
  12. If you use non-standard charges for any residue, change the charges of system 1-3 by changeparm command m on prmtop3, using a pdb file with the correct charges (i.e. the file used as input to pdbtocomqum). Note that leap will change back to standard charges, even if the charges were modified in the input pdb file.

    This is done by command cqcharge, i.e.
    changeparm <<EOF
    prmtop3
    m
    comqum.pdb
    w

    q
    EOF

  13. Normally not needed: Add the desired atom types of the junction atoms into file comqum.dat (if missing)
    after the bond length.
    Sample comqum.dat file.


  14. Unusual: If there is no junctions in the QM system, run command
    cqtrunc_no_junc


  15. Run cqprep and follow the instructions. You will do the following:
    1. Run cqtrunc to generate the file prmtop1 from file prmtop3 using the command
    2. Remove charges from prmtop1 and write out pdbout1 with changeparm (done automatically)
    3. Remove charges of system 1 from prmtop2 and write out pdbout3 with changeparm (done automatically).
    4. Insert a cap into prmtop3 with changeparm - use the data from the equilibration. It is very important that you did not have a cap in prmtop3 before running cqprep.
    5. Run a test of sander.in1 (write down the number of atoms and the energy). Continue the run by writing a q
    6. Run a test of sander.in2 (this may take several minutes; write down the number of atoms and the energy). Continue the run by writing a q

    7. Run checkcoord to see that the correspondance list works (check that there are no warnings (Amber may change the order of the atoms; if so, you have to reorder them also in the original pdb file and rerun everyting from pdbtocomqum (point 3). There should be no warnings. If you get warnings forr HB-atoms, e.g. in CYM, run fixcoord1 (fixcoord1_amberin; fixcoord1_turboin; fixcoord1; fixcoord1_amberout) and see if the program gives errors. If so, you have to swap coordinates in the prmcrd1 file.
    8. Run a test of sander.in3 (write down the number of atoms and the energy

      If you have two junction atoms that are were bound in the full (untruncated system), you need to delete the corresponding bond in the prmtop1 file. This is done with changeparm, command bb on file prmtop1.

  16. Normally not needed: If you run several states and want to compare the energies it is absolutely mandatory that you have exactly the same system2 (i.e. the same belly in sander.in2) in all calculations. You do that in pdbtocomqum by adding or removing residues to system 2 (step 6).
    However, it often enough to copy the same sander.in2 and sander.in3 files to all directories (check with diff that they are identical).


  17. Run define to set up the QM calculation (more instructions; explicit commands):
    1. Make an UFF hessian - coordinate menue, command ff, set the charge and then hit return; when exiting the coordinate menu, you need to explicitly say n, you do not want to define internal coordinates;
    2. Define the basis set.
    3. Set up the wavefunction (eht)
    4. Optionally change the DFT functional
    5. Turn on RI
    6. Turn on dispersion with bj damping.
  18. Remove unneccesary files (an make a backup) by command
    cqzip
    (gzip *; cqgunzip; mkdir Gz; mv *.gz Gz; cp * Gz; gzip Gz/*&)

    The following 15 files are left (for UHF jobs, there are 16 files and alpha and beta replace mos; perhaps also auxbasis):
    basis         control         mos           prmcrd1       prmtop1       prmtop3       sander.in2      sander.out3
    comqum.dat    coord           pointcharges  prmcrd3       prmtop2       sander.in1    sander.in3

  19. Typically you first run a set of protein-fixed calculations (default input; i.e.
    $protein fixed
    is in the file comqum.dat file.
    Then run
    comqum -c 200
    or (with the RI approximation):
    comqum -ri -c 200

    Always check the convergence of your results. We use the standard Turbomole quite floppy convergence criteria, which means that you sometimes can get random convergence. Therefore, check that the energy really is converged.

    The final results are in the prmcrd3 and in the Turbomole files.
    You can build pdb files from those (pdb and pdb3) using the command
    cqtopdb

  20. Then, you run protein-free (i.e. you relax system 2 with MM). Set (in a new directory)
    $protein free
    in file comqum.dat.

    If you have a residue with net charged groups outside the QM system, you need to insert into comqum.dat (change the residue number and the desired charge).
    $residue_charge_corr
     1
     652 -1

    Then, you can start the program with
    comqum -ri -backup -c 200

    Sometimes, the calculation crashes, saying that the charge of the system has changed.
    This normally means that your leap charges (typically of the QM system) were not correct.
    Check what is the correct charge of the protein and if it is not correct in the prmtop3 file, use cqtopdb to write a pdb file with the prmtop charges, then change the charge of one QM atom to get the correct total charge of the protein, and insert those charges into the prmtop3 file with changeparm, command m.

    You can test this part of comqum by running:
    ridft >logd
    sed -i 's/$point_charges/#point_charges/' Str/control
    ridft -proper >mulliken
    sed -i 's/#point_charges/$point_charges/' Str/control
    fixcharge_turboin>>fixc
    fixcharge_amberin>>fixc
    fixcharge>>fixc
    fixcharge_amberout>>fixc

    Sometimes it can be quite a pain to find errors in the charges. Here are some hints.

    Always check the convergence of your results. We use the standard Turbomole quite floppy convergence criteria, which means that you sometimes can get random convergence. Therefore, check that the energy really is converged. This is especially important with protein free.

  21. You should typically run qmmmpbsa on all ComQum results (it takes only ~1 h).
    This is done by
    qmmmpbsa-vac >logv
    qmmmpbsa-ptch>logp
    Note that you need the files, so you cannot copy the calculation to the temp disk of the batch computer.

    This is a first and minimal step to make the QM/MM eneriges more stable (and to include a proper solvation).
    Typically, you should also run calculations forth and back between the start and end states, especially if you have protein free, to ensure that the energies are stable (that the two states reside in the same local minima) or start QM/MM optimisations on several snapshots from a MD simulation.

    Alternatively, you should consider to run QTCP or other QM/MM free-energy perturbation methods.

    If you get the problem
    Too large difference - use $residue_charge_corr

    it often means that you have a different (incorrect) charge of the QM system in the prmtop3 file (for a protein-fixed calculation, this does not matter and is not noticed by comqum until you run qmmmpbsa).
    Add the missing charge to one atom in the QM system in the pdb3 file (by hand) and then run:

    changeparm <<EOF

    prmtop3
    m
    pdb3
    w
    prmtop3
    q
    EOF

    Then, you can run qmmmpdbsa as usual.

  22. You can add restraints on a distance between two atoms using in comqum.dat (note the unusual units):
    $restraints
      atom1(1)  atom2(1)  desired_distance_in_A   force_constant_in_au
    $restraints
     1 3  1.00 1.0
    Note that the atoms should be numbered according to the quantum system and not according to the total (MM) system.
    A typical force constant is 10 a.u. (which typically gives a 0.05 pm deviation at convergence).
    If you get oscillations in the energies or slow convergence, try to reduce the force constant (1.0 a.u. is normally also OK).

    You can also restrain the difference between two distances (r_AB - r_BC):
    $restraints
       atom1 atom2 atom3 desired_distance_in_A   force_constant_in_au
    $restraints
     1 2 3 1.00 1.0


    You can also restrain an angle:
    $restraints
       atom1 atom2 -atom3 desired_angle_in_deg   force_constant_in_au
    Note the minus sign before atom3 (to discern it from the previous restraint).

    $restraints
     1 2 -3 1.00 1.0
    A force constant of 1-10 au seems to be OK.



    You can also restrain a torsion angle:
    $restraints
       atom1 atom2 atom3 atom4 desired_angle_in_deg   force_constant_in_au
    $restraints
     1 2 3 4 1.00 1.0
    Again, a force constant of 1 seems to be OK.

  23. Similarly offset forces are added by (also in comqum.dat):
    $offset forces
       atom1(1) atom2(1) force(1)
       ...
       atom1(n) atom2(n) force(n)
    The atoms should be given by numbers, forces in atomic units (the negative of the output in relax, possibly calculated with
    $interconversion on
       gradient cartesian --> internal
    this works also with redundant internal coordinates if you first define the desired coordinates)
    Note that the atoms should be numbered according to the quantum system and not according to the total system.

  24. In order to calculate strain energies (the first two lines can be done with command str):
    mkdir Str
    cp control coord basis energy mos alpha beta auxbasis Str
    cd Str
    dscf >logd (or ridft >logd)


  25. To evaluate the effect of the protein, you should use the following four data points: The QM/MM energy, E(QM/MM), the QM energy with point charges, E(QM+ptch), found in the fourth column in the energy file, the QM vacuum energy at the QM/MM geometry (optained by the str script or by qmmpbsa-vac, E(QM1), and the QM energy optimised in vacuum, E(QM). Then:
    E(QM1)-E(QM) is the strain energy of the QM system, i.e. the effect of the change in geometry in the protein.
    E(QM+ptch)-E(QM1) is the effect of the point charges in the protein
    E(QM/MM)-E(QM+ptch) is the MM energy. With protein fixed, it is usually small. If not, you should run the reaction forth and back until the energy difference converges. It is typically an artifact of different local minima.

    Typically, the ptch energy dominates, and it can be further understood by using changepdb, command im2, on the two pdb3 files (for two states in an reaction) together with the corresponding charges in Ptch/logd. It will perform a per-residue decomposition of this energy (using a MM approximation, so the difference is not exactly equal to E(QM+ptch)-E(QM1), but normally quite close.

    The MM (=MM12-MM1) energy difference between two states can be understood by running
    changeparm, on parmtop2 command qmd (QM/MM difference)
    Read in the corresponding prmtop1 file (they should apply to both states)
    prmcrd3 and prmcrd1 of the first state
    prmcrd3 and prmcrd1 of the second state

  26. Further improvement of the energy differences can be obtained with big-QM and QTCP calculations.

There is a potential currently unsolved problem in all ComQum calculations:
If junction atoms are involved in improper dihedrals, there will be a mismatch between MM1 and MM3, because the movement of the junction atom will make the improper dihedral angle not equal in MM1 and MM3. It is normally not an important problem. We have seen problems from it only for very distorted structures (caused by other problems).
It is solved by using additive QM/MM.


The ComQum energy is:
    E_QM/MM = E_QM1+ptch23 + E_MM123_no_q1 - E_MM1_no_q
In the energy file, you find the E_QM/MM and E_QM1+ptch23 in the second and fourth columns:
   628 -16161.8463086333  15231.6278196500 -15322.2901877800     -3.4144500185
            E_QM/MM            SCF_kin       E_QM1+ptch23              E_MM_opt
The two MM energies are in the files calcforce.out1 (E_MM1_no_q) and calcforce.out2 (E_MM123_no_q1), which can be accessed by
grep 'Total energy               =' calcforce.out?

Explicit Turbomole define commands

define

y
ff
c -2

no
bb all def2-SV(P)
bl

*
eht

-2
y

dft
func tpss
*
ri
on
m 2500
*
dsp
bj
*
*


How to set up a system with all groups within a specified radius within an original QM system

This can be done with the PMISP tools:
  1. Start from a ComQum calculation of the original QM system (comqum.dat and pdb3 files)
  2. If you have CYM or other strange residues, change them to standard residue names.
  3. Run changepdb on pdb3
  4. Enter command sp
  5. Enter i
  6. Enter the original QM system as ranges (first and last atom; same if no range), with the second atom with a negative number if it is not the last one:
    For example:
    $correspondance_list
        33
       531   534   535   536   537   538   539   540   541  1229  1231  1232
      1233  1234  1266  1269  1270  1271  1272  1273  1274  1275  1276  1324
      1327  1328  1329  1330  1331  1332  1333  1334  1443
    becomes:
     531  
    -531  
     534
    -541
     1229
    -1229
     1231
    -1234
     1266
    -1266
     1269
    -1276
     1324
    -1324
     1327
    -1334
     1443
     1443
  7. Enter the desired radius
  8. Enter return once
  9. Enter auto
  10. Enter restop
  11. Enter return
  12. x
    pdb3.xyz
  13. Finish with q
  14. Finally enter
    pmisp 1 pdb3.xyz split.inp - -
  15. This generates a large number of files. The only you need is
    systa
  16. \rm f*xyz *input prmcrd* comqum_a* c*xyz comqum_b.dat comqum_b1.dat gtot.xyz tot.xyz systb g1.xyz
To avoid so much output, you can also make a file dummy.inp with "0 0" in the first line and an empty second line and then run
pmisp 1 pdb3.xyz split.inp - dummy.inp


If you get when issuing the auto command
Sorry, residue  62 is disconnected.
 This cannot be handled by a standard restop file.
 Please do something about this.
 Quitting...
It typically means that a residues contains two fragments (backbone and sidechain). Try to connect them by adding heavy atoms (in point 8 above).

This is the command cqtrunc
changeparm <<EOF
prmtop3

tr
comqum.dat
n
y
y

prmtop1
y

m
prmcrd3
prmcrd1

q
EOF

In command cqtrunc_no_junct one of the two consecutive y rows is deleted


This is cqprep

# Run cqtrunc
echo "Running cqtrunc"
cqtrunc

# Zero charges in prmtop1
# This is now done by command cqzero1
echo "Zero the charges in prmtop1"
changeparm <<EOF
prmtop1
p
pdbout1
m
prmcrd1
z
w
prmtop1
q
EOF

# Zero the charges of the first QM system in prmtop2
# This is now done by command cqzero3
echo "Zero the charges of syst 1 in prmtop2"
changeparm <<EOF
prmtop3
p
pdbout3
m
prmcrd3
1
w
prmtop2
q
EOF

# Insert the cap into prmtop3
echo " "
echo "Now you should add a cap to prmtop3 if appropriate"
echo "Enter: prmtop3; c; y; enter the appropriate values; w; <CR>; q"
echo " "
changeparm

# Test calcforce
calcforce <<EOF >calcforce.out1
prmtop1
m
prmcrd1
1000
mden1
force1
EOF
less
calcforce.out1
calcforce <<EOF >calcforce.out2
prmtop2
m
prmcrd3
1000
mden2
force2
EOF

less calcforce.out2

# Use define to include the proper basis sets and get a starting wave function.

# DFT functional is already defined
echo "Define basis set and mos by define"
define

# Check the corresponance lists so that Amber has not renumbered the atoms
checkcoord

# Test sander 3 (must be killed)
echo "sander.out3: This job must be killed"
sander -O -i sander.in3 -o sander.out3 -r mdrest3 -e mden3 -c prmcrd3 -p prmtop3
more sander.out3


Sample comqum.dat file

$title
Title
$protein fixed
$junction= 6
$junctions
     1      2     1.09800000000000 H1 
$correspondance_list
15
1976 1978 1979 1980 1981 1982 1983 4865 4866 4867 4868 4869
4870 4871 4872
$junction_atoms
1 2 1.38979963570127
$end

Note that $junctions is added by pdbtocomqum, whereas $junction_atoms is added by changeparm, command tr.

This is a schematic overview of  comqum

# Start with a scf step then loop through the following four steps

# Gradient step (gradient calculation)
grad  >logg
cp gradient grad1
offsetf  >fixf
fixforce_turboin>>fixf
fixforce_amberin>>fixf
fixforce>>fixf
fixforce_turboout>>fixf
mv gradient g1
mv grad1 gradient

# Relaxation step
relax  >logr
fixcoord1_turboin >fix1
fixcoord1_amberin>>fix1
fixcoord1>>fix1
fixcoord1_amberout>>fix1

# Molmech step
# if  protfree then
   sander -O -i sander.in3 -o sander.out3 -r mdrest3 -c prmcrd3 -p prmtop3
   cp mdrest3 prmcrd3

   fixcoord2_turboin >fix2
   fixcoord2_amberin>>fix2
   fixcoord2>>fix2
   fixcoord2_turboout>>fix2
#endif

# Scf step (energy calculation)
calcforce <<EOF >calcforce.out1
prmtop1
m
prmcrd1
1000
mden1
force1
EOF
calcforce <<EOF >calcforce.out2
prmtop2
m
prmcrd3
1000
mden2
force2
EOF

dscf>logd

grep 'Self energy' logd > selfenergy
fixenergy_turboin >fixe
fixenergy_amberin>>fixe
fixenergy>>fixe
fixenergy_turboout>>fixe
#if protfree then
   moloch>mulliken
   fixcharge_turboin >fixc
   fixcharge_amberin>>fixc
   fixcharge>>fixc
   fixcharge_amberout>>fixc
#endif
checkconv



These are the programs you need (apart from Turbomole and Amber):
forcetograd
fixcoord1
fixcoord2
fixenergy
fixcharge
comqum

In addition, you may need some of our accessory programs:
pdbtocomqum
cqprep
espprep
mkchargecorr
checkcoord
cqgunzip
changeparm
changepdb


This is no longer needed (Sept. 2016)
In addition, you need a modified version of sander (sanderf) which writes out the forces.
Three files in Amber/src/sander have been modified (Amber 8.0):
cp md.h md.h_orig
cp force.f force.f_orig
cp mdread.f mdread.f_orig

md.h
8c8
< parameter (BC_MDI=55)
---
> parameter (BC_MDI=54)
28c28
<       idecomp,icnstph,ntcnstph,maxdup,iwforc     !55
---
>       idecomp,icnstph,ntcnstph,maxdup     !54
42c42
<       ivcap,iconstreff,idecomp,klambda,icnstph,ntcnstph,maxdup,iwforc
---
>       ivcap,iconstreff,idecomp,klambda,icnstph,ntcnstph,maxdup

force.f:
five lines added after line 752  (1073 in Amber 9; 550 in Amber 7; 677 in Amber 5; almost at the end);
after the row:
    call trace_exit( 'force' )

   !------- This part added for ComQum ---
   !     Write out the forces
   if (iwforc.eq.1) then
      write(77,'(3e22.14)')(f(i),i=1,3*natom)
      write(77,*)
   endif
   !-----------------------------------

   return
end subroutine force

mdread.f

Three lines as this diff comparison shows:

84c84
<          mmtsb_switch, mmtsb_iterations,rdt,icnstph,solvph,ntcnstph,iwforc &
---
>          mmtsb_switch, mmtsb_iterations,rdt,icnstph,solvph,ntcnstph &
266,268d265 (Before line    !     Check to see if "cntrl" namelist has been defined)
<    !     Write forces?
<           iwforc=0
<




This is an example of a file with the quantum system (syst1)
New title
531
534-541
1229
1231-1234
1266
1269-1276
1324
1327-1334
1443

531
1229
1266
1324

The format is :
A new title line
Atoms of the QM system (one atom or a range on each line)
Blank line
Junction atoms (one on each line)



This is an example of a Gaussian file for calculation of the fitted charges
%Chk=impc42cm.chk
#P HF/Sto-3G SCF(Vshift=1000) Charge=Bohr

Pc.ox/Cu2, ComQum on ImPc42CM, 16/24 A, DZpdf/6-31G*, 3/12-99
pcox.new, Plastocyanin, new equilibration with AMBER, 24 A CAP, 960624

    1    2
h      -2.34543400328041    7.56091572587130   -4.16389870653224
c      -2.51199881855881    8.35999606813360   -3.48399836140879
n      -1.51299928840743    9.17099568670494   -2.99999858904316
c      -2.10299901091925   10.03399528081969   -2.17699897611565
h      -1.58799925313351   10.80799491679282   -1.62799923432075
n      -3.43899838257314    9.80999538617113   -2.10099901185989
h      -4.09199807545487   10.32199514536783   -1.52499928276361
c      -3.72299824900256    8.74799588564985   -2.92699862337644
h      -4.69299779279318    8.29799609729338   -3.07899855188796
h       3.04952150153140    5.77052939610116   -3.07577470019319
c       2.27499893002440    6.42199697961172   -2.68199873860458
h       1.50799929075903    5.80699726885787   -2.21199895965449
h       2.71799872167310    7.06199667860760   -1.91899909745794
s       1.53599927759010    7.43899650129735   -3.92699815305749
h       2.91307908205073   10.89635139119452   -2.44756015205704
c       2.16999897940789   11.41399463177954   -3.03699857164136
n       0.90199957577231   10.89799487446412   -3.20299849356841
c       0.26899987348420   11.77099446387568   -3.94899814271048
h      -0.74299965055302   11.64499452313586   -4.30499797527693
n       1.03399951369021   12.83399396392664   -4.23299800913990
h       0.77999963315122   13.58899360883583   -4.85399771707183
c       2.25999893707918   12.64899405093564   -3.64899828380616
h       3.10999853730807   13.31499373770322   -3.68499826687468
h       2.51873747856117    8.39759526408039    0.78586489906895
c       1.59999924748968    7.82499631975424    0.73999965196398
h       1.82399914213824    6.83299678631063    0.34799983632901
h       1.20799943185471    7.70999637384092    1.75099917647152
s       0.34399983821028    8.58399596278216   -0.25599987959835
c      -0.87399958894124    7.21399660711912   -0.22199989558919
h      -0.42999979776285    6.29999703699063   -0.61599971028353
h      -1.18699944173141    7.02499669600940    0.80499962139325
h      -1.74799917788248    7.47799648295492   -0.81699961574942
cu      0.37499982363039    8.91299580804723   -3.05899856129434

! Point charges
@charges/N

--Link1--
%Chk=impc42cm.chk
#P HF/STO-3G SCF(Vshift=0) Geom=AllCheck Guess=Read Charge=Bohr

! Point charges
@charges/N

--Link1--
%Chk=impc42cm.chk
#P HF/3-21G SCF(Vshift=0) Geom=AllCheck Guess=Read Charge=Bohr

! Point charges
@charges/N

--Link1--
%Chk=impc42cm.chk
#P B3LYP/Gen SCF=(Tight,Vshift=0)  Geom=AllCheck Guess=Read Charge=Bohr

! Gen basis set
@/usr/people/ulf/Basis/stdbas/N

! Point charges
@charges/N

--Link1--
%Chk=impc42cm.chk
#P B3LYP/Gen Geom=AllCheck Guess=(Read,Only) Pop=(Minimal,MK,ReadRadii)
   IOp(6/33=2,6/41=10,6/42=17) Charge=Bohr

! Gen basis set
@/usr/people/ulf/Basis/stdbas/N

! Point charges
@charges/N

! Non-default atomic radii
 Cu 2.0



Reorganisation

Cu2 system

  1. Do the following commands:

  2. mkdir Str Cu1 Cu1p
    cp control coord basis energy mos alpha beta Str
    cp control coord basis energy alpha Cu1
    cp coord basis energy mden1 mden2 prmcrd3 prmtop3 sander.in3 sander.out3 selfenergy pointcharges Cu1p
    cd Cu1
     
  3. Goto Cu1 and in control remove uhf-related keywords and insert
$scfmo    file=mos
$closed shells
 a       1-80                                   ( 2 )
$scfdiis   start=0.5
  1. Also

  2. cp control ../Cu1p
    mv alpha mos
    and insert as the first line in mos
    $scfmo
     
  3. Run the following jobs
cd Str
dscf>logd
cp energy ../Cu1
cd ../Cu1
dscf>logd
cp energy ../Cu1p
cp mos ../Cu1p
cd ../Cu1p
dscf>logd
     
  1. After the jobs have completed run:

  2. cd Str
    mul2
    mv mul ../mul_str
    mv alpha ../alpha_str
    mv beta  ../beta_str
    mv energy ..
    /bin/rm *
    cd ../Cu1
    mul
    mv mul ../mul_cu1
    mv mos ../mos_cu1
    mv energy ..
    /bin/rm *
    cd ../Cu1p
    mul
    fixcharge>fixen
    vi sander.in3
    sander -O -i sander.in3 -o sander.out3 -r mdrest3 -c prmcrd3 -p prmtop3
    fixenergy>>fixen
    mv mul ../mul_cu1p
    mv mos ../mos_cu1p
    mv fixen ..
    mv energy ..
    /bin/rm *
    cd ..
    rmdir Str Cu1 Cu1p




Cu1 system
 
  1. Do the following commands:

  2. mkdir Str Cu2 Cu2p
    cp control coord basis energy mos Str
    cp control coord basis energy mos Cu2
    cp coord basis energy mden1 mden2 prmcrd3 prmtop3 sander.in3 sander.out3 selfenergy pointcharges Cu2p
    cd Cu2
     
     
  3. Goto Cu2 and in control remove uhf-related keywords and insert
$uhf
$uhfmo_alpha    file=alpha
$uhfmo_beta    file=beta
$alpha shells
 a       1-80                                   ( 1 )
$beta shells
 a       1-79                                   ( 1 )
$scfiterlimit       300
$scforbitalshift  closedshell=.05
$natural orbitals    file=natural
$natural orbital occupation    file=natural
  1. Also insert as the first line in mos

  2. $uhfmo_alpha
    Then
    cp control ../Cu2p
    cp mos beta
    mv mos alpha
    and change the first line in beta to
    $uhfmo_beta
     
  3. Run the following jobs
cd Str
dscf>logd
cp energy ../Cu2
cd ../Cu2
dscf>logd
cp energy ../Cu2p
cp alpha ../Cu2p
cp beta ../Cu2p
cd ../Cu2p
dscf>logd
  1. After the jobs have completed run:

  2. cd Str
    mul
    mv mul ../mul_str
    mv mos ../mos_str
    mv energy ..
    /bin/rm *
    cd ../Cu2
    mul2
    mv mul ../mul_cu2
    mv alpha ../alpha_cu2
    mv beta ../beta_cu2
    mv energy ..
    /bin/rm *
    cd ../Cu2p
    mul2
    fixcharge2>fixen
    vi sander.in3
    sander -O -i sander.in3 -o sander.out3 -r mdrest3 -c prmcrd3 -p prmtop3
    fixenergy >>fixen
    mv mul ../mul_cu2p
    mv alpha ../alpha_cu2p
    mv beta ../beta_cu2p
    mv fixen ..
    mv energy ..
    /bin/rm *
    cd ..
    rmdir Str Cu2 Cu2p


 
  • In order to calculate the polarisation energy (Do not do this normally):
    1. Insert into the control file these rows (note that there probably already is a $properties group):
    2. $properties
    3.  potential
    4. $points
    5.  fld
    6.  file = pointcharges
    7. Comment out the first line in file pointcharges:
    8. #$point_charges
    9. Run moloch>fieldfile (it may take quite some time; typically 20 minutes)
    10. Get the polarisation energy by the program polenergy.

    AMBER errors
    Unreasonably large value for MAXPR:
    This indicates that cut=1000 does not work.
    Try to reduce it.


    Old (Amber 5) version:

    How to start a ComQum calculation
      Form to fill in
     
    1. Start with an equilibrated crystal structure (with hydrogen atoms) stored in pdb format with charges.

    2.  
    3. Change the name of metal-bound water molecules to WAM and move them before the other water molecules.

    4.  
    5. Run pdbtocomqum
      1. Select a logfile
      2. Enter the name of the pdb file
      3. Do not search for short contatcts
      4. Enter a new title (can be done in a file, see next point).
      5. Select the quantum system, including the junction atoms that will be changed to hydrogens (you first give the number of the amino acid, then the number of the atoms in this amino acid; one <CR> let you select a new amino acid; another <CR> ends the selection).

      6. This can optionally be done by entering the name of a file with the atom numbers (see below). This file should also contain the new title as the first line.
      7. Select the radius for system 2 (MM system that is optimised). Typically 15 .

      8. If necessary, include or remove residues from this system.
      9. Select the radius for system 3 (MM system that is fixed). Typically the rest of the protein.
      10. Select amber (default)
      11. Select the basis set (i.e. the type of junction, typically 6).
      12. Enter the junction atoms (those that will become hydrogen atoms in the quantum system)

      13. Ensure that they are not unknown.
      14. Remove the charge of neighbours.
      15. Check that the change in the charges is not too large. If so, remove more neighbours.
      16. Dielectric constant of 4 is OK; it is normally not used.

       
    6. Run comqumtoturbo and comqumtoamber

    7. This will give you the following files:
      addfil        control       edit.in1      link.in1      logfil        pdb.in1       pointcharges  sander.in2
      comqum.dat    coord         edit.in3      link.in3      parm.in       pdb.in3       sander.in1    sander.in3
    8. Ensure that all the files are correct.

    9. Normally, only link.in3 needs to be changed:
      1. Include cross-links, if any
      2. Move non-protein residues into a separate molecule for link.in3
      3. Set the protein ends charged if appropriate:
       P   0    0    1    3    1
       
    10. Run cqprep

    11. 1. Check that the number of modified charges in prmtop2 exactly equals the number of atoms in system1.
      2. Check that no electrostatics is present in the run of sander.in1, but in the other two calculations.
      3. Write down the number atoms and the energies in the simulations.
      4. For checkcoord (check if the correspondance list is OK; AMBER may renumber atoms), there should be no warnings (except HB-atoms, e.g. for CYM).
       
    12. If system 1 has a total charge of 0, you should use ESP charges instead of Mulliken charges.
    13. espprep, followed by
      mv tmp.control control
      sets the appropriate keywords in control.
      Otherwise, you should (this technique can also be used on uncharged systems to save time)
      1. Run Gaussian to get an MK ESP fit (or even  RESP ). Note that the point charges must be included in these calculations (copy the file pointcharges to a new name and remove the first and last line; then include this file in the Gaussian calculations).
      2. Run dscf to get the wavefunction
      3. Run moloch to get Mulliken charges (mul2 for a open-shell calculation).
      4. Set $charge_corr by running mkchargecorr

      5.  
    14. It may also be wise to run changeparm, command energy , on prmtop1/3 and prmcrd1/3 too see if there are any problems with the topology.

    15.  
    16. You can add restraints on a distance between two atoms using in comqum.dat (note the unusual units):

    17. $restraints
        atom1(1)  atom2(1)  desired_distance_in_A   force_constant_in_au
      $restraints
       1 3  1.00 10.0
      A typical force constant is 10 a.u. (which typically gives a 0.05 pm deviation at convergence)

      Similarly offset forces are added by (also in comqum.dat):
      $offset forces
         atom1(1) atom2(1) force(1)
         ...
         atom1(n) atom2(n) force(n)
      The atoms should be given by numbers, forces in atomic units (the negative of the output in relax, possibly calculated with
      $interconversion on
         gradient cartesian --> internal
      this works also with redundant internal coordinates if you first define the desired coordinates)

      Note that the atoms should be numbered according to the quantum system and not according to the total system.
       

    18. Remove unneccesary files (an make a backup) by

    19. gzip *; cqgunzip; mkdir Gz; mv *.gz Gz; cp * Gz; gzip Gz/*&
      The following 15 files are left (for UHF jobs, there are 16 files and alpha and beta replace mos)
      basis         control         mos           prmcrd1       prmtop1       prmtop3       sander.in2      sander.out3
      comqum.dat    coord           pointcharges  prmcrd3       prmtop2       sander.in1    sander.in3
       
       
    20. Start comqum with the protein free:

    21. comqum -energy 4 -gcart 2 -c 200
      or (with the RI approximation):
      comqum -ri -energy 4 -gcart 2 -c 200

      It may be necessary to check that fixcharge behaves properly, i.e. that the quantum chemical and original charges are not too different so that the backbone of the junction residues get an erroneous charge.
       
       

    22. After convergence, fix the protein by setting

    23. $protein fixed
      in the control file. Then, restart the job by
      comqum -c 200
      or (with the RI approximation):
      comqum -ri -c 200
      Alternatively, you can run with a fixed protein first (easier to converge and to get a good forceapprox file).
       
    24. In ordet to calculate strain energies:

    25. mkdir Str
      cp control coord basis energy mos alpha beta auxbasis Str
      cd Str
      dscf >logd
      mul or mul2
      mv mul ../mul_str
      mv mos ../mos_str
      mv alpha ../alpha_str
      mv beta ../beta_str
      mv energy ..
      /bin/rm *
      cd ..
      rmdir Str
         

    This is /home/bio/Bin/Comqum/cqprep:
    # Amber system 1
    $OML_AMBER/exe/link -O -i link.in1 -o link.out1 -p $OML_AMBER/dat/db94.dat
    page link.out1
    $OML_AMBER/exe/edit -O -i edit.in1 -o edit.out1 -pi pdb.in1 -po pdbout1
    page edit.out1
    $OML_AMBER/exe/parm -O -i parm.in  -o parm.out1 -f $OML_AMBER/dat/parm96.dat -m $OML_AMBER/dat/Mumod/modparm.dat -p prmtop1 -c prmcrd1
    page parm.out1

    # Amber system 3
    $OML_AMBER/exe/link -O -i link.in3 -o link.out3 -p $OML_AMBER/dat/db94.dat
    page link.out3
    $OML_AMBER/exe/edit -O -i edit.in3 -o edit.out3 -pi pdb.in3 -po pdbout3
    page edit.out3
    $OML_AMBER/exe/parm -O -i parm.in  -o parm.out3 -f $OML_AMBER/dat/parm96.dat -m $OML_AMBER/dat/Mumod/modparm.dat -p prmtop3 -c prmcrd3
    page parm.out3

    # Zero charges in prmtop1 -> prmtop1
    changeparm
    # and in prmtop3 -> prmtop2
    changeparm
    # Add a cap to prmtop3 if appropriate
    changeparm

    # Test sander 1-2
    sanderf -O -i sander.in1 -o sander.out1 -r mdrest1 -e mden1 -c prmcrd1 -p prmtop1
    page sander.out1
    exe/sanderf -O -i sander.in2 -o sander.out2 -r mdrest2 -e mden2 -c prmcrd3 -p prmtop2
    page sander.out2

    # Use define to include the proper basis sets and get a starting wave function.
    # DFT functional is already defined
    define

    # Check the corresponance lists so that Amber has not renumbered the atoms
    checkcoord

    # Test sander 3 (must be killed)
    $OML_AMBER/exe/sander -O -i sander.in3 -o sander.out3 -r mdrest3 -e mden3 -c prmcrd3 -p prmtop3
    page sander.out3


    Sample syst1 file

    The first line is a title
    Then comes all QM atoms as atom numbers or range
    Then a blank line
    Finally the junction atoms, on numer on each line.

    New title
    1210
    1213-1220
    2018
    2020-2023
    2087
    2090-2097
    2145
    2148-2155
    2281

    1210
    2018
    2087
    2145



    Sample leap.in file

    source leaprc.ff99
    loadAmberParams frcmod.ff99SB
    addAtomTypes {{"CX" "C" "sp2"}}
    lloadAmberPrep wam.in
    loadAmberParams wam.dat
    lx=loadpdb pdb.in3
    deleteBond x.174.9 x.175.1
    bond x.771.SG x.796.SG
    saveamberparm x prmtop3 prmcrd3
    quit


    Backup script
    Makes it possible to move files to local disk during execution
    Inserted into comqum the following lines:

    1. after stepgrep() ...}
    # UR Copy files from local to global disk
    backup() {
      if [ "$backup" = "y" ]
      then
          /bin/cp * $PBS_O_WORKDIR
      fi
      }

    2. after     "-ri"    | "ri")    ri="y" ;;
        "-backup" | "backup") backup="y" ;;

    3. After each convgrep at the end of the file (three places)
        checkconv

        backup

    Then run it with the following commands

    mkdir /disk/local/your_user_id
    cd /disk/local/your_user_id
    cp $PBS_O_WORKDIR/* .
    jobex -backup -ri -c 200
    cp * $PBS_O_WORKDIR

    More instructions on hessians

    It is often wise to run a frequency calculation at a lower level of theory (e.g. HF/sto-3g) to get a better forceapprox matrix.
    The easiest thing is to run a uff (MM) calculation with define:
    Go into the molecular geometry menu and select ff
    change the charge of the molecule (if <>0) by c +2
    and then start the calculation with return.
    You have now run a single-point Hessian (frequency) calculation with UFF, and the result is in the file uffhessian0-0.
    Quit define by qq and then change in the control file
    $forceinit on
     diag = carthess
    (instead of default)
    and (re)insert the line
    $grad   file=gradient
    Thereby, you will read the Hessian matrix and use it as an initialisation of the forceapprox matrix.

    Alternatively, if you prefer a higher level of theory (e.g. HF/Sto-3g):
    mkdir Freq
    cp control coord basis energy mos alpha beta auxbasis Freq
    cd Freq

    change basis set to sto-3g hondo and remove dft and ri (with define)

    insert at the end of the control file:
    $dipgrad  file=hessian
    $hessian  file=hessian
    $vibrational normal modes  file=hessian
    $forcepath <path to the scratch files>
    $maxcor <core memory in MB (70% of available memory is recommended >

    dscf >logd
    aoforce >logf  (or if it is an open-shell system: uhfforce >logf)

    hesstofoa
    cp hessfoa ..
    cd ..

    change in the control file
    $forceinit off
    $forceapprox    file=hessfoa

    Obsolete points

    Former point 7: Construct fragment files for system1

    Construct parameter and topology files for the quantum system (without H-atoms).
    These should involve exactly the same atom types and parameters as in the full protein, except at the junctions.


    This is no longer (5/8-06) needed: The parameters are now generated automatically by running changeparm, command tr on the prmtop3 file.
    Go directly to the next point.


    In /home/bio/Data/Comqum/Lubos_database, there is a consistent data base of reasonably truncated variants of all amino acids.

    Old version.
    This can be done automatically by the program mkjunct
    Please, note however, that this program is not thoroughly tested yet and may contain bugs.
    At present, you need to cut out the corresponding parent amino acid from the amber topology file (starting from the name of the amino acid). On the other hand, the program automatically reads the parameter file parm99.dat.
    The output of the program is the new topology and parameter files for the molecule of interest. When tested, these files should be put in /home/bio/Data/Comqum/7 (provided that you use the BP/6-31G* method). These files should be read in into tleap.
    In addition, it gives the file newjunct, which contains the row of the junction. It should be put into the file /home/bio/Data/Comqum/junctfactor in the appropriate place (alphabetical order).
    pdbtocomqum generates the necessary input to mkjunct at the end of the log file (names of junctions and discarded atoms).
    The only additional needed information is the optimum bond length of the junction H atom in the fragment optimised by the same method as in the intended QM/MM calculation.

    Note that the bond parameters of the junction hydrogen atoms are fixed by the corresponding parameters for the junction carbon atoms and junctfact (which is found in the comqum.dat file).

    Thus the force constant of the junction hydrogen atom must be the force constant of the junction carbon atom * juncfact^2, whereas the equilibrium parameter should be the carbon equilibrium parameter / juncfact (i.e. it should be the ideal bond length of that bond with that basis set, as listed in /home/bio/Data/junctfactor_cns. Otherwise a spurious extra force will be introduced:
    FMM1=k(x-x0)^2
    FMM2=k/jf^2(x*jf-x0*jf)^2 = k/jf^2*jf^2(x-x0)^2 = k(x-x0)^2 = FMM1
    Example:
    $junction_atoms
        14    11    1.38765000
    bond CPB  CH3E   572.053     1.5028   (in Amber.dat files)
    bond CPB  HA    1101.53      1.083    (1.083 from QM vacuum opt; 1101.53 calculated on next line:
    1101.53=572.053*(1.38765)^2
    1.38765=1.5028/1.083



    How to prepare for charges for protein free:

    espprep, followed by
    mv tmp.control control
    sets the appropriate keywords in control.

    Otherwise, you should (this technique can also be used on uncharged systems to save time)
    (new procedure by Lubos, according to Markus:)
    Files needed: esp_dscf or esp_ridft, correct_chg
    (if alpha+beta exist you need also: mul2_inp, /home/bio/Bin/mul2)
    1. Make sure that your $TURBODIR points to turbomole version >=5.8
    2. Run script 'esp_ridft' or 'esp_dscf'. It calculates mulliken charges and constructs file comqum.dat.new if everything went well. Copy the section '$charge_corr'  into the  comqum.dat.
    3. change  '$protein fixed' to '$protein free' in comqum.dat
    4. Run'(to get the right charges into prmtop3, and to check that those programs work).
      fixcharge_turboin>>fixc
      fixcharge_amberin>>fixc
      fixcharge>>fixc
      fixcharge_amberout>>fixc
    Older method using only Turbomole calculations (using makepoints and chargefit:
    1. Run dscf/ridft to get the wavefunction
    2. Run makepoints to get proper points in which you will calculate the potential.
    3. Run moloch/moloch2 to get the potentials and also the Mulliken charges (for moloch2, you need to have the file natural and the natural occupation numbers in file control).
    4. Run chargefit to get the ESP charges (there are many choices here; use the default but restrain to the potential with the weight 1).
    5. Run mkchargecorr to get the $charge_corr

      Yong Shen has set up automatic shells for this. They can be found in /home/bio/Sheny/script/protein_free1/protein_free
      or on toto7/whenim64
      /home/sheny/script/protein_free(1)
      They need to be slightly modified for your personal use.

      Old procedure:
    1. Run Gaussian to get an MK ESP fit (or even  RESP ). Note that the point charges must be included in these calculations (copy the file pointcharges to a new name and remove the first and last line; then include this file in the Gaussian calculations).
    2. Run dscf to get the wavefunction
    3. Run moloch to get Mulliken charges (mul2 for a open-shell calculation).
    4. Set $charge_corr by running  mkchargecorr

    Running without cut=1000

    Recent results indicates that you can get more stable convergence if you change cut=1000 in sander.in3 (if you have enough memory for it). Then you can run comqum with normal convergence (as in the previous point) and need not to fix the protein again.

    Otherwise:
    Run
    comqum -energy 4 -gcart 2 -c 200
    or (with the RI approximation):
    comqum -ri -energy 4 -gcart 2 -c 200

    It may be necessary to check that fixcharge behaves properly, i.e. that the quantum chemical and original charges are not too different so that the backbone of the junction residues get an erroneous charge.

    After convergence, fix the protein by setting

    $protein fixed
    in the comqum.dat file. Then, restart the job by
    comqum -c 200
    or (with the RI approximation):
    comqum -ri -c 200
    Alternatively, you can run with a fixed protein first (easier to converge and to get a good forceapprox file).

    With the (obsolete) version 2, it is often seen that the energy increases abruptly when protein is set free. This comes from the self-energies which are incorrectly calculated by the QM calculations (1-2 and 1-3 interactions are included).

    With version 3, this should no happen.


    How to start an old calculation

    1. cp control comqum.dat
    2. vi comqum.dat
      A. delete non-comqum rows (keep $title, $restraingts, $protein fixed/free, $junction, $junction atoms, $correspondance_list)
      B. Add at the end
      $version
       02
      $end
      C. Also add the number of atoms as the first row after $correspondance_list
    3. Delete the following keywords in all three sander.in files:
      init=3
      nrun=1
      cut2nd=20.0
      idiel=1

      Also change cut=1000, in sander.in1 and sander.in2


    Script to automatically run many constrained simulation (Hao 11/3-21)


    #!/bin/sh

    #SBATCH -t 144:00:00

    #SBATCH -N 1

    #SBATCH -n 32

    #SBATCH --exclusive

    #SBATCH -A SNIC2020-3-18

     

    export PARA_ARCH=SMP

    export PARNODES=32

    export PATH=$TURBODIR/scripts:$TURBODIR/bin/x86_64-unknown-linux-gnu_smp:$PATH

     

    A=162

    B=176

    initdist=2.00

    finaldist=1.00

    steplength=-0.1

    bond=Fe2-H

    #

    sed -i '/\$end/i\$restraints' comqum.dat;sed -i "/\$end/i\ $A $B $initdist 1.0" comqum.dat

    cd $SNIC_TMP

    cp -p $SLURM_SUBMIT_DIR/* .

    comqum -backup -c 600 -ri

    cp -pu * $SLURM_SUBMIT_DIR

    cd $SLURM_SUBMIT_DIR

    Purtur

    mkdir $bond-$initdist

    cp * ./$bond-$initdist

    #

    while [ $(echo "$initdist != $finaldist" | bc) -eq 1 ]

    do

    initdist=`echo $initdist + $steplength |bc`

    sed -i '$d' comqum.dat

    sed -i '$d' comqum.dat

    sed -i "/\$restraints/a\ $A $B $initdist 1.0" comqum.dat

    sed -i '$a\$end' comqum.dat

    cd $SNIC_TMP

    cp -p $SLURM_SUBMIT_DIR/* .

    comqum -backup -c 600 -ri

    cp -pu * $SLURM_SUBMIT_DIR

    cd $SLURM_SUBMIT_DIR

    Purtur

    mkdir $bond-$initdist

    cp * ./$bond-$initdist

    done


    History

    5/6-03: Changes in the self energy of systems 2-3
    Before: self energy calculated by QM program and added to total energy
    After:   self energy calculated by MM2 program
    Why: otherwise 1,2 and 1,3 electrostatic interactions are included in self energy. Therefore, the energy differs from what you get from a MM calculation. Moreover, the energy contains informatin from system 3 that should not change.
    Note that we still keep an infinite cutoff for the interaction between system 1 and the other subsystems.

    Specific changes
    1. fixenergy_turboin.f: self energy removed
    2. fixenergy_turboout.f: self energy removed
    3. cqprep: command 1 (zero system 1 charges), instead of z (zero all charges) in prmtop2
    4. changeparm: new command 1 (zero system 1 charges, using $correspondance in comqum.dat)
    5. No: changed back again (16/6): comqumtoamber.f: nsnb=1,cut=25 in sander.in1 and sander.in2; no belly in sander.in2.
    6. 13/10-03: sander.in1/2 cut=1000 (otherwise, the energies become unstable and strange) (comqumtoamber)
    7. 13/10-03: both energies can be calculated by ComQum: in comqum.dat $version 02 or 03 (default) (fixenergy, fixenergy.turboin, fixenergy.turboout)
    How to recalculate the energy to go from version 2 to version 3
    1. Change in comqum.dat
      $version
      03

    2. Set cut=1000 in sander.in1 and in sander.in2.

    3. Run changeparm on prmtop 3, command 1 (zero charges of syst 1, the program then reads comqum.dat) and then w(rite) to prmtop2. Save copies of the old files. Enter:

      changeparm
      prmtop3
      1
      w
      prmtop2
      q

    4. Run
      sander -O -i sander.in1 -o sander.out1 -r mdrest1 -e mden1 -c prmcrd1 -p prmtop1
      \rm forcedump.dat
      sander -O -i sander.in2 -o sander.out2 -r mdrest2 -e mden2 -c prmcrd3 -p prmtop2
      \rm forcedump.dat


    5. Replace the last row in file energy with that from file job.last (almost at the end of the file):

      Last energy change:
      Old:     947.366120 H
      New:       0.000000 H

         126 -4320.99279891900  4303.70087538700 -8624.69367430600                                   Copy this row
         126 -5268.35892311563  4303.70087538700 -5269.53383013889    -9.13165371658

      If job.last is deleted, you can recalculate this row in the following way:
      The first energy is the third ComQum energy plus the self energy (run calcselfen to get it).
      The second energy is also the second ComQum energy
      The third energy is the first energy minus the second energy.

      126 cq3-selfen cq2 cq3-selfen-cq2

      126 cq1        cq2 cq3             cq4

      In the example above selfen=-948.5410312199, so
      -5269.53383013889 -  -948.5410312199 =
      -4320.99279891900
       4303.70087538700                    =  4303.70087538700
      -4320.99279891900 - 4303.70087538700 = -8624.69367430600

      Alternatively, you can simply rerun dscf (ridft) again.

    6. Run
      fixenergy_amberin
         >fixen
      fixenergy_turboin  >>fixen

      fixenergy          >>fixen
      fixenergy_turboout >>fixen

    7. Now you should have a new corrected total energy as the final row in file energy and a detailed account of the energies in file fixen.
    How to recalculate the energy to go from version 3 to version 2
    1. Change in comqum.dat
      $version
      02

    2. Run changeparm

      changeparm
      prmtop3
      z
      w
      prmtop2
      q

    3.  Set cut=15 in sander.in1 and in sander.in2.

    4. Run
      sander -O -i sander.in1 -o sander.out1 -r mdrest1 -e mden1 -c prmcrd1 -p prmtop1
      \rm forcedump.dat
      sander -O -i sander.in2 -o sander.out2 -r mdrest2 -e mden2 -c prmcrd3 -p prmtop2
      \rm forcedump.dat

    5. Replace the last row in file energy with that from file job.last (almost at the end of the file):

      Last energy change:
      Old:     947.366120 H
      New:       0.000000 H

         126 -4320.99279891900  4303.70087538700 -8624.69367430600                                   Copy this row
         126 -5268.35892311563  4303.70087538700 -5269.53383013889    -9.13165371658

      If job.last is deleted, you can recalculate this row in the following way:
      The first energy is the third ComQum energy.
      The second energy is also the second ComQum energy
      The third energy is the first energy minus the second energy.

      126 cq3  cq2 cq3-cq2

      126 cq1  cq2 cq3      cq4

      Alternatively, you can simply rerun dscf (ridft) again.

    6. Run
      fixenergy_amberin
         >fixen
      fixenergy_turboin  >>fixen

      fixenergy          >>fixen
      fixenergy_turboout >>fixen

    7. Now you should have a new corrected total energy as the final row in file energy and a detailed account of the energies in file fixen.
    If you get problem to run sander.in2 (too large pairlist)
    1. Set in comqum.dat
      $version
      02

    2. Run changeparm on prmtop3, zero charges, and write it back to prmtop2
      changeparm
      prmtop3
      z
      w
      prmtop2
      q

    3. Set cut=15 in sander.in1 and in sander.in2.
    Then energies will be calculated the old way with self energy taken from QM calculation. Normally, the difference is small.

    040224: Fixed error in fixcharge.f
    Added the line (48):
          rfirst(nresi+1)=natom3
    Errors only if the last residue is in the quantum system and protein is free.

    041110: Fixed error in fixenergy.f
    Moved (line)57
    call CqRdVersion(file,version,.true.)
    outside the
    if (found_elstat) then
    section
    Before this correction, version 2 energies do not contain the selfenergy.

    160527: Inserted cqtrunc into cqprep

    160907: Retired sanderf; use instead dumpfrc

    sanderf -O -i sander.in1 -o sander.out1 -r mdrest1 -e mden1 -c prmcrd1 -p prmtop1
    mv fort.77 force1
    sanderf -O -i sander.in2 -o sander.out2 -r mdrest2 -e mden2 -c prmcrd3 -p prmtop2
    mv fort.77 force2

    Old sander.in1
    sander.in1                                         
     &cntrl
      irest=0,ntx=1,
      nstlim=1,dt=0.0,
      imin=0,maxcyc=1,drms=0.001,
      temp0=300.0,ntt=1,tautp=0.2,
      ntc=1,ntf=1,
      nsnb=25,cut=1000.0,dielc=1.0,
      ntpr=1,ntwx=0,ntwv=0,ntwe=1,
      ntb=0,ntp=0,taup=0.2,
      iwforc=1,ibelly=0
     &end

    Old sander.in2
    Comqum syst1 file auto-generated by pmisp for mph                              
     &cntrl
      irest=0,ntx=1,
      nstlim=1,dt=0.0,
      imin=0,maxcyc=1,drms=0.001,
      temp0=300.0,ntt=1,tautp=0.2,
      ntc=1,ntf=1,
      nsnb=25,cut=1000.0,dielc=1.0,
      ntpr=1,ntwx=0,ntwv=0,ntwe=1,
      ntb=0,ntp=0,taup=0.2,
      iwforc=1,ibelly=0
     &end

    14/11-16: Introduced calcforce (instead of sander)
    Retired dumpfrc (it gave erroneous results) and sander

    sander -O -i sander.in1 -o sander.out1 -r mdrest1 -e mden1 -c prmcrd1 -p prmtop1
    more sander.out1
    \mv forcedump.dat force1
    \cp force1 mden1
    sander -O -i sander.in2 -o sander.out2 -r mdrest2 -e mden2 -c prmcrd3 -p prmtop2
    more sander.out2
    \mv forcedump.dat force2

    \cp force2 mden2

    MoFe cluster in N2ase, Neutralised protein, 2H state, subunit C, LC & UR 5/9-16
     &cntrl
      irest=0,ntx=1,
      nstlim=1,dt=0.0,
      imin=0,maxcyc=1,drms=0.001,
      temp0=300.0,ntt=1,tautp=0.2,
      ntc=1,ntf=1,
      nsnb=25,cut=1000.0,dielc=1.0,
      ntpr=1,ntwx=0,ntwv=0,ntwe=1,ntxo=1,
      ntb=0,ntp=0,taup=0.2,
      ibelly=0
     &end

     &debugf
      do_debugf=1,dumpfrc=1
     &end

    MoFe cluster in N2ase, Neutralised protein, 2H state, subunit C, LC & UR 5/9-16
     &cntrl
      irest=0,ntx=1,
      nstlim=1,dt=0.0,
      imin=0,maxcyc=1,drms=0.001,
      temp0=300.0,ntt=1,tautp=0.2,
      ntc=1,ntf=1,
      nsnb=25,cut=30.0,dielc=1.0,
      ntpr=1,ntwx=0,ntwv=0,ntwe=1,ntxo=1,
      ntb=0,ntp=0,taup=0.2,
      ibelly=0
     &end

     &debugf
      do_debugf=1,dumpfrc=1
     &end

    3/1-17 Implemented additive and automatic subtractive QM/MM,

    as well as the possibility to ignore a constant MM3 term with protein fixed.

    They are started by the following keywords in comqum.dat
    $subtractive - runs a subtractive QM/MM, using the files prmtop1/2, prmcrd1/3, mden1/3, force1/3
    $additive - runs an additive QM/MM, using the files prmtop2, prmcrd3, mden1/3, force1/3 (mden1 and force1 are dummy files with zeroed energies and forces)
    $fixed mm3 - tells the program to calculate the MM-only terms only once. They will be reported in the keyword
    $mm3_constant
       57154
      0.149838102904E+04  0.885893514365E+04  0.458467360620E+05  0.181372676476E+05
      0.164620408024E+06  0.102008364009E+06 -0.115430580572E+07 -0.813335713801E+06
    The first number is the number of atoms. The others are the bond, angle, dihedral, 1,4 van der Waals, 1,4 electrostatics, other van der Waals, other electrostatics, and total energy, respectively in kJ/mol.

    Hopefully, this will be transparent to the users (i.e. the program without extra keywords in comqum.dat will run as before).

    To run with $additive, you need to set up the prmtop1 and prmcrd1 (with cqtrunc), as well as dummy mden1 and force1 files with zero results before starting the calculation:

    mden1:
    L0  Nsteps           time(ps)         Etot             EKinetic
    L1  Temp             T_solute         T_solv           Pres_scal_solu
    L2  Pres_scal_solv   BoxX             BoxY             BoxZ
    L3  volume           pres_X           pres_Y           pres_Z
    L4  Pressure         EKCoM_x          EKCoM_y          EKCoM_z
    L5  EKComTot         VIRIAL_x         VIRIAL_y         VIRIAL_z
    L6  VIRIAL_tot       E_pot            E_vdw            E_el
    L7  E_hbon           E_bon            E_angle          E_dih
    L8  E_14vdw          E_14el           E_const          E_pol
    L9  AV_permMoment    AV_indMoment     AV_totMoment     Density          dV/dlambda
    L0        1  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
    L1  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
    L2  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
    L3  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
    L4  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
    L5  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
    L6  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
    L7  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
    L8  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00
    L9  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00  0.0000000000E+00

    force1:
    Run program dumforce and enter the number of atoms
    or copy this line natom times:
      0.00000000000000E+00  0.00000000000000E+00  0.00000000000000E+00

    May 2021: Parallel version of calcforce (by Pedro Ojeda May).
    calcforce-omp.f
    So far implemented only in comqum-par

    On Kebnekaise
    export TURBODIR=/proj/nobackup/teobio/Bio/TURBO/Turbo7.5.1
    export PATH=$TURBODIR/scripts:$TURBODIR/bin/em64t-unknown-linux-gnu_mpi:$PATH
    export PARA_ARCH=SMP
    export PARNODES=28
    export PATH=$PATH:/proj/nobackup/teobio/Bio/Bin

    module purge
    module add GCC/8.2.0-2.31.1
    module load foss/2019a
    ulimit -s unlimited
    export OMP_NUM_THREADS=8

    ./comqum-par -c 200 -ri

    On Tetralith
    Compile with
    module load buildtool-easybuild/4.3.3-nscf4a9479 foss/2019a
    \rm ../Utilib/*.o
    make calcforce-omp
    \rm ../Utilib/*.o
    But do not load this module when you run:

    ulimit -s unlimited
    export OMP_NUM_THREADS=8
    ./comqum-par -c 200 -ri




    Tidy up after comqum in many directories (if purtur is not run):
    gzip -r *
    for x in GEO_OPT_CONVERGED.gz GEO_OPT_RUNNING.gz fixcharge.out.gz fixcoord1.out.gz fixcoord2.out.gz fixenergy.out.gz fixforce.out.gz job.start.gz mden.gz mden1.gz mden2.gz mdinfo.gz mdrest1.gz mdrest2.gz mdrest3.gz mulliken.gz nextstep.gz not.converged.gz pbs_save_files.gz pulf.gz restartg.gz sander.out1.gz sander.out2.gz statistics.gz temp.ulf.gz truc.gz uffhessian0-0.gz ; do
     find * -name $x -delete -print
    done