I'm getting a weird problem when trying to minimize a solvated protein-ligand complex, in which I treat the ligand as a QM region. Well, the program stops since at some point SCCDFTB cannot anymore converge, but the key problem happens a bit before. My ligand is a zwitter-ion, having a COO- and NH3+ groups. The NH3 group is by default oriented towards a nearby ASP residue. And at the very first steps of minimization the ASP shifts closer to the ligand's NH3 group so that there is a covalent-bond-like distance between ASP's O and ligand's N atoms, and then the O atom totally overlaps with one of the NH3 group's H atoms. It seems like an oxygen atom doesn't "feel" the qm protons (or another way around, since electrostatics is treated on qm part), and it also looks as if this NH3 group proton is considered to be a linking atom, but there are no QM/MM bonded interactions in my definition, therefore no linking atoms...
I use the following commands to initiate SCCDFTB: SCCDFTB remove CHRG 0 - sele segid ACPC end - TEMP 100.00 SCFT 0.0000001 - DISP MULL
And it seems that Charmm understands it quite well:
CHARMM> SCCDFTB remove CHRG 0 - CHARMM> sele segid ACPC end - CHARMM> TEMP 100.00 SCFT 0.0000001 - CHARMM> DISP MULL SELRPN> 14 atoms have been selected out of 24264 SCCINT> REMOve: Classical energies within QM atoms are removed. SCCINT> No EXGRoup: QM/MM Elec. for link atom host only is removed. SCCINT> No QINP: Charges will be based on atomic numbers. Dispersion among QM atoms included SCCINT> DO NOT USE CUT-OFF FOR QM/MM Mulliken population requested ------------------------------------------------ SCCINT: Classical atoms excluded from the QM calculation: NONE. SCCINT: Quantum mechanical atoms: 4589 ACPC 1 ACPC C2 4590 ACPC 1 ACPC H21 4591 ACPC 1 ACPC H22 4592 ACPC 1 ACPC C3 4593 ACPC 1 ACPC H31 4594 ACPC 1 ACPC H32 4595 ACPC 1 ACPC C1 4596 ACPC 1 ACPC N 4597 ACPC 1 ACPC HT1 4598 ACPC 1 ACPC HT2 4599 ACPC 1 ACPC HT3 4600 ACPC 1 ACPC C4 4601 ACPC 1 ACPC OT1 4602 ACPC 1 ACPC OT2 SCCINT: Quantum mechanical link atoms: NONE. SCCINT: Quantum mechanical GHO boundary atoms: NONE. ------------------------------------------------ FINDEL: Quantum atom 4589 ACPC 1 ACPC C2 assigned to element: C 6 FINDEL: Quantum atom 4590 ACPC 1 ACPC H21 assigned to element: H 1 FINDEL: Quantum atom 4591 ACPC 1 ACPC H22 assigned to element: H 1 FINDEL: Quantum atom 4592 ACPC 1 ACPC C3 assigned to element: C 6 FINDEL: Quantum atom 4593 ACPC 1 ACPC H31 assigned to element: H 1 FINDEL: Quantum atom 4594 ACPC 1 ACPC H32 assigned to element: H 1 FINDEL: Quantum atom 4595 ACPC 1 ACPC C1 assigned to element: C 6 FINDEL: Quantum atom 4596 ACPC 1 ACPC N assigned to element: N 7 FINDEL: Quantum atom 4597 ACPC 1 ACPC HT1 assigned to element: H 1 FINDEL: Quantum atom 4598 ACPC 1 ACPC HT2 assigned to element: H 1 FINDEL: Quantum atom 4599 ACPC 1 ACPC HT3 assigned to element: H 1 FINDEL: Quantum atom 4600 ACPC 1 ACPC C4 assigned to element: C 6 FINDEL: Quantum atom 4601 ACPC 1 ACPC OT1 assigned to element: O 8 FINDEL: Quantum atom 4602 ACPC 1 ACPC OT2 assigned to element: O 8
SCCDFN> Some atoms will be treated quantum mechanically.
The number of SCCDFTB QM atoms = 14 The number of molecular mechanical atoms = 24250 The number of MM atoms excluded from QM = 0 Of which the number of QM/MM link atoms = 0
so I am confused, why should such an overlap happen. The whole set of input/output files are in the attachment (the log is in separate archive, since it is too big).
I will greatly appreciate any comments/hints of what can be the reason of the problem.
In fact, another thing worries me in the log file. After initial declaration of nonbonded interactions parameters Charmm calculates reasonable number of Atom pairs :
General atom nonbond list generation found: 12060632 ATOM PAIRS WERE FOUND FOR ATOM LIST 702420 GROUP PAIRS REQUIRED ATOM SEARCHES
Then, after SCCDFTB parameters were set and "update" command issued, it prints:
There are 12060632 atom pairs and 44839 atom exclusions. There are 0 group pairs and 4104 group exclusions. Generating nonbond list with Exclusion mode = 5 == PRIMARY == SPACE FOR 15687725 ATOM PAIRS AND 0 GROUP PAIRS NBONDS: Using routine NBONDA for list generation.
General atom nonbond list generation found: 40376 ATOM PAIRS WERE FOUND FOR ATOM LIST 2286 GROUP PAIRS REQUIRED ATOM SEARCHES
Well, may be here Charmm has recalculated only atom pairs for the QM region (although I doubt so, since the system size is more then 20000 atoms, and 40376 number seems to be too small assuming, there is no cut-off for QM/MM nonbonded interactions), but then, starting minimization (which as we already know doesn't lead to anything good), it prints :
NONBOND OPTION FLAGS: ELEC VDW ATOMs CDIElec FSHIft VATOm VFSWIt BYGRoup NOEXtnd NOEWald CUTNB = 14.000 CTEXNB =999.000 CTONNB = 8.000 CTOFNB = 12.000 WMIN = 1.500 WRNMXD = 0.500 E14FAC = 1.000 EPS = 1.000 NBXMOD = 5 There are 40376 atom pairs and 44839 atom exclusions. There are 0 group pairs and 4104 group exclusions. Generating nonbond list with Exclusion mode = 5 == PRIMARY == SPACE FOR 15687725 ATOM PAIRS AND 0 GROUP PAIRS NBONDS: Using routine NBONDA for list generation.
General atom nonbond list generation found: 40376 ATOM PAIRS WERE FOUND FOR ATOM LIST 2286 GROUP PAIRS REQUIRED ATOM SEARCHES
But where are all those another millions of atom pairs?
As allways, thanks in advance for your kind comments)
OK, I solved the problem. It seems, that although I have manually replaced "1" column for "0" column in psf file for the atoms that were kept frozen in the preliminary solvent equilibration, this was not enough. Typing "cons fix sele none end" at the beginning of the script solved the problem.
