$molecule
3 1
Al 0.00000000000000 0.00000000000000 0.00000000000000
O 1.80185960291498 -0.00000982662585 -0.22421140724593
H 2.51639802970329 0.87356529080086 -0.12246105621416
H 2.51638889484236 -0.87359208188566 -0.12245800808720
$end
$rem
JOBTYPE FREQ
METHOD HF
VIBMAN_PRINT 6
NO_REORIENT TRUE
MAX_SCF_CYCLES 100
BASIS Aug-cc-pvtz
SCF_CONVERGENCE 10
THRESH 14
$end
But the results that I obtained were exactly the same. So I am not sure if the MBD contribution is being considerate. Do you know if there is a way of putting the mbd contribution on the second derivative?
The bug that you’ve discovered is that MBD frequencies aren’t supported. The code should have told you this, and defaulted to finite difference, but it didn’t. I will submit a ticket.
I checked, and both IDERIV=0 and IDERIV=1 give the same result. There are two issues, which are that (1) the code should automatically go to finite-diff, and (2) apparently the MBD energy and gradient info aren’t added to the quantities that are differenced in the finite-diff code.
My question is, if I want to see than the MBD contribution to the force constant, is it correct if I just subtract the HF-MBD force constant - HF force constant?
Those differences (sub-wavenumber) are in the convergence noise, I think. I don’t believe it’s possible to compute MBD frequencies with any current version of Q-Chem.
Hi John, So I was doing some tests and tried to obtain the gradients itself for the MBD (I ran an HF-MBD and PBE0-MBD) with this two inputs:
$molecule
3 1
Al 0.00500000000000 0.00000000000000 0.00000000000000
O 1.80185960291498 -0.00000982662585 -0.22421140724593
H 2.51639802970329 0.87356529080086 -0.12246105621416
H 2.51638889484236 -0.87359208188566 -0.12245800808720
$end
$rem
JOBTYPE force
METHOD HF
NO_REORIENT TRUE
MAX_SCF_CYCLES 100
BASIS aug-cc-pVTZ
MBDVDW MBD_FORCES
MBDVDW_BETA 10000
SCF_CONVERGENCE 10
THRESH 14
$end
$molecule
3 1
Al 0.00500000000000 0.00000000000000 0.00000000000000
O 1.80185960291498 -0.00000982662585 -0.22421140724593
H 2.51639802970329 0.87356529080086 -0.12246105621416
H 2.51638889484236 -0.87359208188566 -0.12245800808720
$end
$rem
JOBTYPE force
METHOD PBE0
NO_REORIENT TRUE
MAX_SCF_CYCLES 100
BASIS aug-cc-pVTZ
MBDVDW MBD_FORCES
MBDVDW_BETA 8500
SCF_CONVERGENCE 10
THRESH 14
$end
But for both outputs I got the same error on the freq calculation for the MBD gradients:
Energy from LibMBD: -6.767440482398079E-004
Force calculation requested from LibMBD, resulting forces in a.u.
Atom 1 : NaN NaN
NaN
Atom 2 : NaN NaN
NaN
Atom 3 : NaN NaN
NaN
Atom 4 : NaN NaN
NaN
But I used to calculate the gradients for MBD all the time with Qchem 6.2.2 (see output ex below):
Force calculation requested from LibMBD, resulting forces in a.u.
Atom 1 : 5.131402789079820E-006 1.342466005685693E-005
-3.325264617820776E-005
Atom 2 : -1.138175572790299E-005 2.207709259368315E-005
-2.893676216140745E-005
Atom 3 : 1.926679125458893E-005 1.492121304558350E-005
-3.344003762037749E-005
Atom 4 : -1.135398386864196E-005 -8.084185932355496E-006
-4.021526120510593E-005
Atom 5 : -8.099570615636543E-007 -2.960170487256917E-005
9.394213031909857E-005
Atom 6 : -1.780055515351012E-006 2.218217141953931E-005
-7.213803477315047E-005
Atom 7 : -2.436613481808033E-005 -1.664334197878195E-005
3.710773661039494E-005
Atom 8 : 1.728280549705957E-005 -2.929746094132065E-005
3.293562194378651E-005
Atom 9 : 8.010887450811649E-006 1.102155660936424E-005
4.399725306496917E-005
+------------------------------------------+
| END Calculate MBD vdW contributions |
+------------------------------------------+
Do you think it might be something with the new version? or my input is wrong?
No, there was a bug introduced in the MBD gradient at some point, it was noted also here:
and I opened a ticket already. I think something got broken in work with Tkatchenko’s group to add new libMBD features to Q-Chem but I’m not sure exactly.
I see, but does this bug appeared recently? because I was calculating 1 week ago with 6.2 the gradients. So I don’t know when it appeared or if its broken also on the older versions.
Hi John, I tested the MBD forces on both methods (with and without finite-diff) with the old version of qchem (6.2) and they gave me similar results. Than I tested again the FREQ calculation with finite-diff to see if something changed now that the MBD forces are computed. So I compared again the single HF with the HF-MBD with these two inputs:
$molecule
3 1
Al 0.00000000000000 0.00000000000000 0.00000000000000
O 1.80185960291498 -0.00000982662585 -0.22421140724593
H 2.51639802970329 0.87356529080086 -0.12246105621416
H 2.51638889484236 -0.87359208188566 -0.12245800808720
$end
$rem
JOBTYPE FREQ
IDERIV 0
METHOD HF
VIBMAN_PRINT 6
NO_REORIENT TRUE
MAX_SCF_CYCLES 100
BASIS Aug-cc-pvtz
SCF_CONVERGENCE 10
THRESH 14
$end
$molecule
3 1
Al 0.00000000000000 0.00000000000000 0.00000000000000
O 1.80185960291498 -0.00000982662585 -0.22421140724593
H 2.51639802970329 0.87356529080086 -0.12246105621416
H 2.51638889484236 -0.87359208188566 -0.12245800808720
$end
$rem
JOBTYPE FORCE
METHOD HF
NO_REORIENT TRUE
MBDVDW MBD_FORCES
MBDVDW_BETA 10000
MAX_SCF_CYCLES 100
BASIS Aug-cc-pvtz
SCF_CONVERGENCE 10
THRESH 14
$end
So this time I got some larger difference on the output:
Only HF:
Vibrational Frequencies in atomic units
-0.0000000000 -0.0000000000 0.0000000000 0.0000000000 0.0000000000 0.0000000000
0.0143452775 0.0484443606 0.0498777267 0.1114262341 0.1333656444 0.1474101892
Vibrational Electric Transition dipole: 1 0.00000000000000 0.00000000000000 0.00000000000000
Vibrational Electric Transition dipole: 2 0.00000000000000 0.00000000000000 0.00000000000000
Vibrational Electric Transition dipole: 3 0.00000000000000 0.00000000000000 0.00000000000000
Vibrational Electric Transition dipole: 4 0.00000000000000 0.00000000000000 0.00000000000000
Vibrational Electric Transition dipole: 5 0.00000000000000 0.00000000000000 0.00000000000000
Vibrational Electric Transition dipole: 6 0.00000000000000 0.00000000000000 0.00000000000000
**********************************************************************
** **
** VIBRATIONAL ANALYSIS **
** -------------------- **
** **
** VIBRATIONAL FREQUENCIES (CM**-1) AND NORMAL MODES **
** FORCE CONSTANTS (mDYN/ANGSTROM) AND REDUCED MASSES (AMU) **
** **
**********************************************************************
Mode: 1 2 3
Frequency: 615.69 1131.43 1148.05
Force Cnst: 2.3597 0.9367 0.8742
Red. Mass: 10.5655 1.2420 1.1257
IR Active: YES YES YES
Raman Active: YES YES YES
X Y Z X Y Z X Y Z
Al -0.392 0.000 0.048 0.011 -0.000 0.020 -0.000 -0.030 -0.000
O 0.600 -0.000 -0.089 -0.010 0.000 -0.121 0.000 0.079 0.000
H 0.484 -0.012 0.058 -0.065 -0.029 0.698 0.656 -0.222 -0.130
H 0.484 0.012 0.058 -0.065 0.029 0.698 -0.656 -0.222 0.130
TransDip 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
Mode: 4 5 6
Frequency: 1715.93 1877.27 1973.65
Force Cnst: 1.8856 2.3453 2.3448
Red. Mass: 1.0869 1.1295 1.0217
IR Active: YES YES YES
Raman Active: YES YES YES
X Y Z X Y Z X Y Z
Al -0.005 0.000 -0.001 0.000 0.007 -0.000 0.002 0.000 -0.000
O -0.072 0.000 -0.003 -0.000 -0.090 0.000 0.030 -0.000 0.006
H 0.641 -0.292 0.037 0.319 0.620 0.096 -0.268 -0.653 -0.041
H 0.641 0.292 0.037 -0.319 0.620 -0.096 -0.268 0.653 -0.041
TransDip 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
STANDARD THERMODYNAMIC QUANTITIES AT 298.15 K AND 1.00 ATM
HF-MBD:
brational Frequencies in atomic units
-0.0095308677 -0.0000000000 -0.0000000000 -0.0000000000 0.0000000000 0.0000000000
0.0000000000 0.0186908969 0.0395645335 0.1480780125 0.1967981655 0.2851517078
Vibrational Electric Transition dipole: 1 0.00000000000000 0.00000000000000 0.00000000000000
Vibrational Electric Transition dipole: 2 0.00000000000000 0.00000000000000 0.00000000000000
Vibrational Electric Transition dipole: 3 0.00000000000000 0.00000000000000 0.00000000000000
Vibrational Electric Transition dipole: 4 0.00000000000000 0.00000000000000 0.00000000000000
Vibrational Electric Transition dipole: 5 0.00000000000000 0.00000000000000 0.00000000000000
Vibrational Electric Transition dipole: 6 0.00000000000000 0.00000000000000 0.00000000000000
**********************************************************************
** **
** VIBRATIONAL ANALYSIS **
** -------------------- **
** **
** VIBRATIONAL FREQUENCIES (CM**-1) AND NORMAL MODES **
** FORCE CONSTANTS (mDYN/ANGSTROM) AND REDUCED MASSES (AMU) **
** **
**********************************************************************
Mode: 1 2 3
Frequency: -501.85 702.78 1022.49
Force Cnst: 0.1769 4.4866 0.7318
Red. Mass: 1.1923 15.4180 1.1880
IR Active: YES YES YES
Raman Active: YES YES YES
X Y Z X Y Z X Y Z
Al -0.016 -0.021 -0.013 -0.456 -0.007 0.062 -0.025 0.024 -0.010
O -0.002 0.066 0.081 0.760 0.016 -0.127 0.016 -0.072 0.066
H 0.619 -0.303 -0.568 0.246 0.138 0.152 -0.240 0.153 -0.288
H -0.160 -0.172 -0.368 -0.109 -0.191 0.207 0.668 0.344 -0.510
TransDip 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
Mode: 4 5 6
Frequency: 1978.11 2280.43 2745.01
Force Cnst: 2.3635 3.3007 4.8087
Red. Mass: 1.0252 1.0773 1.0832
IR Active: YES YES YES
Raman Active: YES YES YES
X Y Z X Y Z X Y Z
Al -0.008 -0.001 0.005 -0.004 0.001 0.001 -0.019 -0.001 0.007
O -0.005 -0.014 -0.028 -0.024 -0.063 -0.008 -0.041 0.031 -0.041
H 0.278 -0.525 0.166 0.614 0.735 0.113 0.399 -0.120 0.203
H 0.020 0.773 0.144 -0.113 0.231 -0.007 0.762 -0.363 0.263
TransDip 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
STANDARD THERMODYNAMIC QUANTITIES AT 298.15 K AND 1.00 ATM
Do you think now it can being computed the MBD contribution to the total force constants?
I’ve looked into this some more, and I now think the vdW-TS energy does get added into the total SCF energy, so finite-diff with IDERIV=0 should give correct frequencies. If those agree with IDERIV=1 (which it seems to do in my tests), then those are also correct frequencies. I’ve noticed that the finite-diff procedure is a bit unusually sensitive to thresholds and I suspect that has to do with computing the Hirshfeld volumes that are needed for the dispersion correction. The analytic 2nd derivatives are definitely not implemented and that should get backstopped but does not. (I will add for next release.)
Note: I made this a new topic because the discussion has strayed pretty far from your original question about DC-DFT.