QM/MM calculations

Aarzoo · September 27, 2024, 5:35am

Hi
I want to do QM/MM calculation. A total of 36 molecules and 1 center molecule is the QM part and the remaining 35 molecules are the MM part. I want to fix the MM part, only the QM part will relax.
This is my input:
$molecule
0 1
coordinates
$end

$rem
basis = cc-pVDZ
METHOD = wB97XD
qm_mm_interface = janus
force_field = amber99
user_connect = true
jobtype = opt
symmetry = false
symmetry_ignore = true
thresh = 14
SCF_CONVERGENCE = 8
$end

$force_field_params
NumAtomTypes 6
AtomType -1 -0.02950000 3.2058099 0.2042
AtomType -2 0.05170000 3.3152123 0.0988
AtomType -3 0.00260000 3.3152123 0.0988
AtomType -4 -0.07780000 3.3152123 0.0988
AtomType -5 0.13950000 2.6254785 0.0161
AtomType -6 -0.07550000 3.3152123 0.0988
Bond -1 -2 610.60 1.4210
Bond -2 -3 904.32 1.3510
Bond -2 -4 573.02 1.4760
Bond -4 -5 791.44 1.0860
Bond -4 -6 736.88 1.4060
Angle -2 -1 -4 129.38 126.610
Angle -4 -1 -4 134.52 120.050
Angle -3 -2 -1 174.43 119.190
Angle -3 -2 -4 131.47 127.520
Angle -4 -2 -1 174.43 119.190
Angle -4 -3 -4 131.47 127.520
Angle -4 -3 -2 131.48 127.520
Angle -2 -4 -4 133.16 120.820
Angle -4 -4 -4 137.53 120.020
Angle -4 -4 -5 97.36 119.880
Angle -4 -6 -6 132.24 121.110
Angle -4 -6 -4 137.75 118.380
Angle -6 -6 -4 132.24 121.110
Torsion -1 -4 -6 -6 3.625 180.0 2
Torsion -1 -2 -3 -4 6.650 180.0 2
Torsion -1 -2 -4 -4 0.700 180.0 2
Torsion -1 -4 -6 -4 3.625 180.0 2
Torsion -1 -4 -4 -4 3.625 180.0 2
Torsion -2 -1 -4 -6 0.300 180.0 2
Torsion -2 -1 -4 -4 0.300 180.0 2
Torsion -2 -3 -4 -4 0.505 180.0 2
Torsion -2 -4 -4 -4 3.625 180.0 2
Torsion -4 -1 -2 -3 0.728 180.0 2
Torsion -3 -2 -4 -4 0.505 180.0 2
Torsion -4 -2 -3 -4 8.510 180.0 2
Torsion -4 -4 -4 -4 3.625 180.0 2
Torsion -4 -3 -4 -4 0.505 180.0 2
Torsion -3 -4 -4 -4 3.625 180.0 2
Torsion -4 -1 -2 -4 0.728 180.0 2
Torsion -4 -4 -3 -4 0.505 180.0 2
Torsion -4 -1 -4 -6 0.300 180.0 2
Torsion -4 -6 -4 -4 3.625 180.0 2
Torsion -4 -6 -6 -4 0.795 180.0 2
Torsion -4 -1 -4 -4 0.300 180.0 2
Torsion -4 -4 -6 -4 3.625 180.0 2
Torsion -4 -4 -6 -6 3.625 180.0 2
Torsion -4 -4 -4 -6 3.625 180.0 2
Torsion -6 -6 -4 -4 3.625 180.0 2
Torsion -2 -1 -4 -4 1.100 180.0 2
Torsion -1 -2 -3 -4 1.100 180.0 2
Torsion -4 -4 -3 -2 1.100 180.0 2
Torsion -4 -4 -4 -2 1.100 180.0 2
Torsion -4 -4 -4 -3 1.100 180.0 2
Torsion -1 -4 -6 -4 1.100 180.0 2
Torsion -4 -4 -6 -6 1.100 180.0 2
Torsion -1 -4 -4 -5 3.625 180.0 2
Torsion -2 -4 -4 -5 3.625 180.0 2
Torsion -3 -4 -4 -5 3.625 180.0 2
Torsion -4 -4 -4 -5 3.625 180.0 2
Torsion -5 -4 -4 -4 3.625 180.0 2
Torsion -5 -4 -4 -5 3.625 180.0 2
Torsion -4 -6 -4 -5 3.625 180.0 2
Torsion -4 -4 -4 -5 1.100 180.0 2
Torsion -4 -6 -4 -5 1.100 180.0 2
Torsion -1 -4 -4 -6 1.100 180.0 2
$end

$qm_atoms
1401:1456
$end

$opt
FIXED
1:1400 XYZ
1457:2016 XYZ
ENDFIXED
$end

But In output:

                   Maximum     Tolerance    Cnvgd?
     Gradient           0.030641      0.000300      NO
     Displacement       0.068384      0.001200      NO
     Energy change           NaN      0.000001      NO

Energy is NaN

In the output, the energy value is not coming. what is this NaN?

kaushik · September 27, 2024, 6:47am

NaN stands for “not a number”, perhaps resulting from two or more atoms or charges being too close to each other such that the corresponding interaction energy blows up. You can perhaps double check your coordinates again. Also run an MM-only calculation to isolate the cause. If these checks are okay, add “NO_REORIENT = true” to your QM/MM calculation’s REM. If this doesn’t fix the issue, post the geometry here or email it to Q-Chem Support.

jherbert · September 27, 2024, 1:12pm

I agree with Kaushik’s suggestion but also I really think you’d be better off running these calculations through the Amber/Q-Chem, GROMACS/Q-Chem, or CHARMM/Q-Chem interfaces.

Aarzoo · September 30, 2024, 12:55pm

Hi @kaushik and @jherbert
Thanks for the suggestions. There was some mistake in defining the atom connectivity and now it is resolved and the calculation was completed at the DFT level.
I want to do spin-flip TDDFT. For QM/MM SF-TDDFT I made the input:

$molecule
0 3
coordinates
$end

$rem
basis = cc-pVDZ
METHOD = wB97XD
qm_mm_interface = janus
force_field = amber99
user_connect = true
jobtype = opt
symmetry = false
symmetry_ignore = true
thresh = 14
MAX_CIS_CYCLES = 700
SCF_CONVERGENCE = 8
MAX_SCF_CYCLES = 700
geom_opt_MAX_CYCLES = 900
UNRESTRICTED = true
model_system_mult = 3
sts_mom = true
SPIN_FLIP = 1
CIS_N_ROOTS = 4
CIS_STATE_DERIVATIVE = 1
$end

$force_field_params
NumAtomTypes 6
AtomType -1 -0.02950000 3.2058099 0.2042
AtomType -2 0.05170000 3.3152123 0.0988
AtomType -3 0.00260000 3.3152123 0.0988
AtomType -4 -0.07780000 3.3152123 0.0988
AtomType -5 0.13950000 2.6254785 0.0161
AtomType -6 -0.07550000 3.3152123 0.0988
Bond -1 -2 610.60 1.4210
Bond -2 -3 904.32 1.3510
Bond -2 -4 573.02 1.4760
Bond -4 -5 791.44 1.0860
Bond -4 -6 736.88 1.4060
Angle -2 -1 -4 129.38 126.610
Angle -4 -1 -4 134.52 120.050
Angle -3 -2 -1 174.43 119.190
Angle -3 -2 -4 131.47 127.520
Angle -4 -2 -1 174.43 119.190
Angle -4 -3 -4 131.47 127.520
Angle -4 -3 -2 131.48 127.520
Angle -2 -4 -4 133.16 120.820
Angle -4 -4 -4 137.53 120.020
Angle -4 -4 -5 97.36 119.880
Angle -4 -6 -6 132.24 121.110
Angle -4 -6 -4 137.75 118.380
Angle -6 -6 -4 132.24 121.110
Torsion -1 -4 -6 -6 3.625 180.0 2
Torsion -1 -2 -3 -4 6.650 180.0 2
Torsion -1 -2 -4 -4 0.700 180.0 2
Torsion -1 -4 -6 -4 3.625 180.0 2
Torsion -1 -4 -4 -4 3.625 180.0 2
Torsion -2 -1 -4 -6 0.300 180.0 2
Torsion -2 -1 -4 -4 0.300 180.0 2
Torsion -2 -3 -4 -4 0.505 180.0 2
Torsion -2 -4 -4 -4 3.625 180.0 2
Torsion -4 -1 -2 -3 0.728 180.0 2
Torsion -3 -2 -4 -4 0.505 180.0 2
Torsion -4 -2 -3 -4 8.510 180.0 2
Torsion -4 -4 -4 -4 3.625 180.0 2
Torsion -4 -3 -4 -4 0.505 180.0 2
Torsion -3 -4 -4 -4 3.625 180.0 2
Torsion -4 -1 -2 -4 0.728 180.0 2
Torsion -4 -4 -3 -4 0.505 180.0 2
Torsion -4 -1 -4 -6 0.300 180.0 2
Torsion -4 -6 -4 -4 3.625 180.0 2
Torsion -4 -6 -6 -4 0.795 180.0 2
Torsion -4 -1 -4 -4 0.300 180.0 2
Torsion -4 -4 -6 -4 3.625 180.0 2
Torsion -4 -4 -6 -6 3.625 180.0 2
Torsion -4 -4 -4 -6 3.625 180.0 2
Torsion -6 -6 -4 -4 3.625 180.0 2
Torsion -2 -1 -4 -4 1.100 180.0 2
Torsion -1 -2 -3 -4 1.100 180.0 2
Torsion -4 -4 -3 -2 1.100 180.0 2
Torsion -4 -4 -4 -2 1.100 180.0 2
Torsion -4 -4 -4 -3 1.100 180.0 2
Torsion -1 -4 -6 -4 1.100 180.0 2
Torsion -4 -4 -6 -6 1.100 180.0 2
Torsion -1 -4 -4 -5 3.625 180.0 2
Torsion -2 -4 -4 -5 3.625 180.0 2
Torsion -3 -4 -4 -5 3.625 180.0 2
Torsion -4 -4 -4 -5 3.625 180.0 2
Torsion -5 -4 -4 -4 3.625 180.0 2
Torsion -5 -4 -4 -5 3.625 180.0 2
Torsion -4 -6 -4 -5 3.625 180.0 2
Torsion -4 -4 -4 -5 1.100 180.0 2
Torsion -4 -6 -4 -5 1.100 180.0 2
Torsion -1 -4 -4 -6 1.100 180.0 2
$end

$qm_atoms
1401:1456
$end

$opt
FIXED
1:1400 XYZ
1457:2016 XYZ
ENDFIXED
$end

Here I have given the overall multiplicity 3 and model_system_mult = 3.

If I am not defining this model_system_mult = 3 in $rem then it gives a multiplicity error. So, is this input is correct?

jherbert · September 30, 2024, 1:21pm

What are you trying to accomplish? For singlet excited states using spin flip, then a triplet is the correct reference state.

Aarzoo · September 30, 2024, 1:36pm

I want to optimize my geometry at E.S=1 using SF-TDDFT. Now I corrected cis_state_derivative = 1 in input.

jherbert · September 30, 2024, 3:50pm

I don’t understand what “E.S=1” means. For a triplet reference state, the singlet ground state (S0) will be the first excited state. Therefore, if your intention is to optimize on S1 (first singlet excited state), then you want CIS_STATE_DERIV=2.

Aarzoo · October 1, 2024, 4:24am

I intend to optimize the S0 state (cis_state_deriv=1).