Comqum

Setup of ComQum calculations
The comqum program is described in

A complete is also given in http://signe.teokem.lu.se/~ulf/Methods/comqum.html

Prepare syst1 file with QM system
A syst1 file contains the QM system. It is structures as

Title

xx--yy pdb atom numbers for atoms in the QM region.

x junction atoms

y

An example is given below

1-27 # First residue in QM region (numbers refer to pdb file) 

1259 # First atom of next residue in QM region (junction atom) 

1261-1271 # here comes the residue 

2576 # more residues (junction)  2578-2592 # more... </tt>

3358-3364 # This is a metal and a few waters </tt>

27 # Junctions - Caps are always CH3 groups junctions should be between C-C bonds </tt>

1259 </tt>

2576 </tt>

Run pdbtocomqum

 * 1) type pdbtocomqum in a terminal
 * 2) select logfile
 * 3) n (no short contacts)
 * 4) enter title
 * 5) Specify QM system (best done by file "syst1", see below for the syntax of this file)
 * 6) Enter cut off radius for system 2 (normally 6 is fine). System 2 is the MM system that is optimized
 * 7) N (do not remove amino acids from system 2)
 * 8) N (do not include more amino acids in system 2)
 * 9) Enter cut-off radius for system 3 (1000 is fine). System 3 is the MM system that is fixed
 * 10) use new format
 * 11) Enter junction atoms. This can also be through the syst1 file. Check that the junctions are known (i.e. not -1)
 * 12) Do not remove charge (n)
 * 13) Do not redistribute charges (n)

Making turbomole and amber files
Run comqumtoturbo</tt> and comqumtoamber</tt>

Make partop3 prmcrd3 with leap
copy the leap.in file (and all .in and .dat files + the leapprc file) from the equilibrium run and modify it:

x=loadpdb pdb.in3</tt>

saveamberparm x prmtop3 prmcrd3</tt>

Delete the line with solventcap (e.g. like solvatecap x TIP3PBOX {0.0 0.0 0.0}  40).</tt>

Then run

tleap -s -f leap.in</tt>

A sample file is give below (make sample file at the bottom)

Typical problems in leap
Leap connects everything that is not separated by "TER". Therefore, make sure to separate groups that should not be connected with TER in the input file

Modify non-standard residues
If there was non-standard residues in the equilibrium calculations, they should be modified, also in the comQum calculation. Use the script

changeparm <<EOF</tt>

prmtop3</tt>

m</tt>

comqum.pdb</tt>

w</tt>

<tt>q</tt>

<tt>EOF</tt>

Add junctions to comqum.dat
continue from here

Run cqprep
Type cqprep and follow the instructions (a cap is inserted - use the data from the equilibrium run).

1) prmtop1 is generated from prmtop3. What is also done: In prmtop1 is a (rudementary) force field for the QM system, so that the MM energy of the QM system can be calculated (given in sander.out1).

2) Charges from prmtop1 are removed and pdbout1 is written out (what is in this file?)

3) Charges from system 1 are removed from prmtop2 and pdbout3 is written (what is in this file?)

4) A cap is inserted (is done manually - follow the instructions)

5) Test run sander.in1 (write down energies)

6) Test run sander.in2 (write down energies)

7) Define is started (select DFT method, bases etc. manually)

Clean-up
gzip *; cqgunzip; mkdir Gz; mv *.gz Gz; cp * Gz; gzip Gz/*&

Run with protein free
1) Replace <tt>$protein fixed</tt> with <tt>$protein free</tt> in comqum.dat

2) Insert <tt>$esp_fit kollman</tt> in the control file

3) Start comqum normally

4) In case of problems run:

<tt>dscf >logd (or ridft >logd)</tt>

<tt>ridft - proper mulliken (comment out the pointcharge line in the control file)</tt>

<tt>fixcharge58_turboin >> fixc</tt>

<tt>fixcharge_amberin >> fixc</tt>

<tt>fixcharge>>fixc</tt>

5) The are sometimes problems with the charges (one problem can occur if the QM system is slightly different from the system used to obtained the RESP charges for non-standard residues). This can be solved by <tt>fixcharge_amberin</tt>
 * For all residues where there are differences between the current QM system and the system used to obtained charges, replace the charges of these residues with standard AMBER charges (see e.g. another residue of the same amino acid in the pdb file) the system pdb file
 * Insert the new charges in the prmtop3 file with changeparm
 * Run

<tt>fixcharge_turboin</tt>

<tt>fixcharge</tt>

<tt>fixcharge_amberout</tt>

Note: Comqum expects that all residues are neutral (as done in the AMBER force field). Some force fields woks with non-neutral residues (e.g. the GLYCAM force field). If some residues are not neutral, it can be necessary to add a correction via comqum.dat. This should only be done with care! normally, a error in the fixcharge is rather due to a problem in the setup. An example of where a charge-correction is necessary is when a force field, e.g., for a substrate is defined as several residues with different charge. In this case, a charge correction can be necessary if the QM region only includes some of the residues in the substrate (since ComQum expects netural residues).

<tt>$residue_charge_corr</tt>

<tt>n # numer of residue to be corrected</tt>

<tt>m c # m= residue id, c = correction </tt>

For example

<tt>$residue_charge_corr</tt>

<tt>1</tt>

<tt>240 0.1940</tt>

Run then

<tt>fixcharge_amberin</tt>

<tt>fixcharge_turboin</tt>

<tt>fixcharge</tt>

<tt>fixcharge_amberout</tt> ￼ == MM analysis == In some (rare) cases it is possible to get a large differences by energies in the MM parts: This is often (but not always) caused by vdW terms. Changeparm can be used to analyse the contributions from individual terms. This is done as follows:

Collect the comqum calculations in two directories called dir1 and dir2

<tt>changeparm</tt>

<tt>path_to_dir1/prmtop2</tt>

<tt>qmd # This is the analysis command</tt>

<tt>path_to_dir1/prmtop1</tt>

<tt>m</tt>

<tt>path_to_dir1/prmcrd3</tt>

<tt>m</tt>

<tt>path_to_dir1/prmcrd1</tt>

<tt>m</tt>

<tt>path_to_dir2/prmcrd3</tt>

<tt>m</tt>

<tt>path_to_dir2/prmcrd1</tt>

Differences in vdW are expected to be around 0-1 kc/mol - much larger values calls fror investigation.