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...
3358-3364 # This is a metal and a few waters
27 # Junctions - Caps are always CH3 groups junctions should be between C-C bonds
1259
2576
Run pdbtocomqum[]
- type pdbtocomqum in a terminal
- select logfile
- n (no short contacts)
- enter title
- Specify QM system (best done by file "syst1", see below for the syntax of this file)
- Enter cut off radius for system 2 (normally 6 is fine). System 2 is the MM system that is optimized
- N (do not remove amino acids from system 2)
- N (do not include more amino acids in system 2)
- Enter cut-off radius for system 3 (1000 is fine). System 3 is the MM system that is fixed
- use new format
- Enter junction atoms. This can also be through the syst1 file. Check that the junctions are known (i.e. not -1)
- Do not remove charge (n)
- Do not redistribute charges (n)
Making turbomole and amber files[]
Run comqumtoturbo and comqumtoamber
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
saveamberparm x prmtop3 prmcrd3
Delete the line with solventcap (e.g. like solvatecap x TIP3PBOX {0.0 0.0 0.0} 40).
Then run
tleap -s -f leap.in
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
prmtop3
m
comqum.pdb
w
q
EOF
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 $protein fixed with $protein free in comqum.dat
2) Insert $esp_fit kollman in the control file
3) Start comqum normally
4) In case of problems run:
dscf >logd (or ridft >logd)
ridft - proper mulliken (comment out the pointcharge line in the control file)
fixcharge58_turboin >> fixc
fixcharge_amberin >> fixc
fixcharge>>fixc
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
- 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
fixcharge_amberin
fixcharge_turboin
fixcharge
fixcharge_amberout
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).
$residue_charge_corr
n # numer of residue to be corrected
m c # m= residue id, c = correction
For example
$residue_charge_corr
1
240 0.1940
Run then
fixcharge_amberin
fixcharge_turboin
fixcharge
fixcharge_amberout

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
changeparm
path_to_dir1/prmtop2
qmd # This is the analysis command
path_to_dir1/prmtop1
m
path_to_dir1/prmcrd3
m
path_to_dir1/prmcrd1
m
path_to_dir2/prmcrd3
m
path_to_dir2/prmcrd1
Differences in vdW are expected to be around 0-1 kc/mol - much larger values calls fror investigation.