Result analysis of QM/MM output from Gromacs_CP2K

QM/MM for Biomolecular Simulation
42

Hello Dmitri,

I am trying to run QM/MM simulation according to your tutorial. But I am facing some charge distribution problems. Gromacs is showing me a warning:-
##########################################################################
QMMM Note: topologly in .tpr will be modified
Number of MM atoms= 6800; Number of QM atoms= 79
Total MM charge= 0.30580
Charge removed from QM atoms= -0.30580
Warning: Total charge of you QMMM system is 0.00000 (QM) + 0.30580 (MM) = 0.30580,
which is not the same as it was in classical MD system!
In order to keep total charge of system in QMMM simulations,
consider to manualy spread -0.30580 charge over MM atoms nearby to the QM region

######################################################################

Kindly help to resolve the charge problem.

Thanking you,
Ranabir

43

yes this is normal situation, it just tells you that total QMMM system charge will differ a little from MM. There is no any particular particular way to resolve it. One of that is suggested in the message from GROMACS. On practice however, I never saw problems with that excess charge.

44

Dear Dmitry,

Thanks a lot for your kind advice and clearing my doubts.

With regards,’
Ranabir

45

Dear Dmitry,

I want to tell you that I successfully ran the TDDFT calculation. But I am now slightly confused about the HOMO-LUMO gap and TDDFT excitation analysis. I ran my QM/MM calculation with the gmx_cp2k package for 100 steps. I have got the two HOMO-LUMO gaps and five excitation energy states in the TDDFT excitation analysis part for every step in cp2k.out file.

Kindly suggest to me, which HOMO-LUMO gap and excitation state I should consider and from which step of the QM/MM simulation.

For the purpose mentioned above, I am attaching below a part from the cp2k.out file. Kindly take a look.

With regards,
Ranabir

HOMO - LUMO gap [eV] : 0.270678
HOMO - LUMO gap [eV] : 1.567375

** ######## ####### ####### ######## ####### ######## **
** ## ## ## ## ## ## ## ## ## **
** ## ## ## ## ## ###### ####### ## **
** ## ## ## ## ## ## ## ## **
** ## ####### ####### ## ## ## **

  •                        TDDFPT Initial Guess                             -

      State         Occupied      ->      Virtual          Excitation
      number         orbital              orbital          energy (eV)

         1              124     (alp)        125 (alp)        0.27068
         2              123     (bet)        124 (bet)        1.56737
         3              123     (alp)        125 (alp)        1.56807
         4              124     (alp)        126 (alp)        1.75016
         5              123     (bet)        125 (bet)        1.93415

  Number of active states:                                     159191

  •                  TDDFPT WAVEFUNCTION OPTIMIZATION                       -

      Step            Time         Convergence           Conv. states

         1            49.8          1.1175E-02                      2
         2            40.9          1.1148E-02                      3
         3            44.9          1.0858E-03                      4
         4            40.7          2.4549E-04                      4
         5            36.4          2.3629E-05                      5

------------------------- Restart Davidson iterations -------------------------
6 34.2 3.3706E-15 5

  • TDDFPT run converged in 6 iteration(s)

U-TDDFPT states of multiplicity 2

     State    Excitation        Transition dipole (a.u.)        Oscillator
     number   energy (eV)       x           y           z     strength (a.u.)
     ------------------------------------------------------------------------

TDDFPT| 1 0.27084 5.1283E-02 4.0716E-02 2.2530E-02 3.18197E-05
TDDFPT| 2 1.56705 3.8467E-02 7.5957E-02 7.0921E-02 4.71412E-04
TDDFPT| 3 1.56915 4.1398E-02 4.7431E-02 6.0582E-02 2.93463E-04
TDDFPT| 4 1.86429 3.4183E-01 7.1103E-01 1.6012E-01 2.95993E-02
TDDFPT| 5 1.93470 2.6487E-02 2.0252E-02 5.3633E-02 1.89037E-04

TDDFPT : CheckSum = 0.128412E+00

  •                        Excitation analysis                              -

    State             Occupied              Virtual             Excitation
    number             orbital              orbital             amplitude

         1   0.27084 eV
                           124 (alp)            125 (alp)        -0.999997
         2   1.56705 eV
                           123 (bet)            124 (bet)         0.815224
                           123 (alp)            125 (alp)        -0.579125
         3   1.56915 eV
                           123 (alp)            125 (alp)         0.815206
                           123 (bet)            124 (bet)         0.579100
         4   1.86429 eV
                           124 (alp)            126 (alp)        -0.980140
                           121 (bet)            126 (bet)        -0.065179
                           119 (bet)            126 (bet)         0.063173
                           124 (alp)            130 (alp)        -0.062887
         5   1.93470 eV
                           123 (bet)            125 (bet)         0.999656

Decoupling Energy: 0.0000379555
Recoupling Energy: -0.0001007892

ENERGY| Total FORCE_EVAL ( QMMM ) energy [a.u.]: -423.510205278879994

46

Could you upload somewhere both cp2k input an output files?

47

Dear Dmitry,
Kindly check the input, log, and out files in the below link:-

https://drive.google.com/drive/folders/11jAlUQzaoxGqvdZU8JXqEBm6MduroYhI?usp=sharing

Thanking you,
Ranabir Majumder

48

Hi Ranabir,

How does system looks like? Are you sure that you have upaired electron in that?

Best,
Dmitry

49

Dear Dmitry,

I have checked the system in gaussian view. There, I have found its charge (0) and multiplicity (2).

50

Could you upload md-qm.pdb file ?

51

Dear Dmitry,
Kindly check the md-qm.pdb file below.

md-qm.pdb (325.2 KB)

52

This is not possible, by default all such kind of biological molecules should have multiplicity (1). Probably your ligand has some charge, check again system. What is this ligand molecule? could you give me a name?

53

Dear Dmitry,

It is a complex system of an organic molecule (Riboflavin) and a DNA base (G).

54

It should not have non-paired electrons by default, so something is wrong with you total QM charge.