i am doing TDDFT calculations with the PCM module (IEFPCM to be precise) and was wondering if it is possible to place User-specified cavities around or in the vicinity of my molecule? In a sense, such a “dummy atom” would then cut out a portion of the continuum but not alter the TDDFT calculation itself as if a real atom would be placed there. I’ve been looking at the $solvent / “SolventAtom” keyword but didn’t manage to come up with a clear solution.
I would be very happy for any comments or suggestions.
You can certainly introduce user-defined atomic radii for a vdw-type cavity. I tried to do this with a ghost atom (thinking you could then use a user-defined ghost atom radius), but it appears that even if I try to assign a radius for atomic number Z=0 (ghost atom) in the $van_der_waals section, that the PCM code zeroes out that radius. Therefore I think the answer is no.
I have one idea on a possible workaround with the PEqS code; will update here if it works.
Thank you very much for your reply! I have just tried your suggestion to introduce a user-defined radius for a ghost atom Z=0 and it seems to modify the amount of tesserae and molecular surface area correctly in my minimal example. I only get issued the warning that “No atomic polarizability available.” which seems clear for the ghost atoms. I’ve uploaded my 2 In-/Output-Files for easier comparison here (one calculation with and without ghost atoms).
If you say that the PCM code zeroes out the radius, i suppose the difference is then due to change in basis set when “ghost atoms” are introduced?
Yes, I’m pretty sure that it’s not working as you intend it. Set PCM_PRINT = 2 in $rem to get more info about the cavity construction, e.g., for glycine with a ghost atom
SWITCHING FUNCTION PARAMETERS (Angst)
Atom R(sphere) R(inner) R(switch)
$molecule
0 1
C -3.036207 0.829905 0.115962
C -1.509446 0.929890 0.362587
O -1.069144 2.079579 0.650443
O -0.874422 -0.148361 0.269649
N -3.621214 2.212792 0.084392
H -3.506594 0.278525 0.932459
H -3.255316 0.328086 -0.826731
H -2.896823 2.832271 0.493661
H -4.491605 2.300730 0.615716
H -3.802442 2.538915 -0.869226 @F 0.0 0.0 0.0
$end
Furthermore, with PCM_PRINT = 2 you get the cavity discretization grid points output in PQR format, which you can copy/paste out of Q-Chem output file into .pqr file that can be read by VMD. If you do that, it’s quite clear that the ghost atom is not contributing to the cavity surface.
I also checked whether this could be done with the Poisson Equation Solver (PEqS) code (https://doi.org/10.1063/1.5023916) but it seems that it has the same issue. In hindsight this is not surprising because PEqS uses the same subroutines to construct a cavity. With the Q-Chem 5.4.1 bugfix release there will be a way to input arbitrary cavity surface (defined pointwise on a grid), although it’s kind of a pain - we did it for a collaborator who wanted multiple dielectric regions (https://doi.org/10.1021/acs.jpcb.0c04032).
If you really want to do this, send me an email and maybe we can collaborate on something. I think it’s a relatively easy fix (famous last words…). However, it occurs to me that one may need to be a little bit careful to avoid possible artifacts of TDDFT. Won’t know for sure without trying it.
First of all, thank you for putting so much effort into my question! I just ran your example molecule (glycine with a ghost atom) on my machine with QChem 5.0.1 to check out the PCM_PRINT=2 statement and get a different output than you
I’ve also checked the cavity discretization grid points output in PQR format as you suggested and with VMD it seems that the ghost atom does indeed contribute to the cavity surface (see Fig. below). Is it possible that there was a change from version 5.0.1 to 5.4? Am I completely off here and doing something completely wrong?
I’ve uploaded both in- & output again here for easier comparison.
$comment
PCM on glycine with a ghost atom
$end
$molecule
0 1
C -3.036207 0.829905 0.115962
C -1.509446 0.929890 0.362587
O -1.069144 2.079579 0.650443
O -0.874422 -0.148361 0.269649
N -3.621214 2.212792 0.084392
H -3.506594 0.278525 0.932459
H -3.255316 0.328086 -0.826731
H -2.896823 2.832271 0.493661
H -4.491605 2.300730 0.615716
H -3.802442 2.538915 -0.869226
@F 0.0 0.0 0.0
$end
$rem
sym_ignore true
EXCHANGE camb3lyp
BASIS 6-31G*
SOLVENT_METHOD pcm
CIS_N_ROOTS 1
CIS_SINGLETS 1
CIS_TRIPLETS 0
PCM_PRINT 2
$end
$pcm
theory iefpcm
method swig
printlevel 2
radii read
$end
$solvent
dielectric 2.27
$end
$van_der_waals
1
0 3.0
$end
and the complete output-file I get is uploaded in the link in my previous reply. I also ran the computation with QChem 4.4, which we also have on one of our local machines and it seems that the cavities for ghost atoms can be assigned there as well (although I had to specify the cavity radius in the $van_der_waals section also for each non-ghost atom, respectively, which was not necessary in QChem 5.0.1).
Anyways, thanks a lot for going through this with me, your suggestions really helped me a ton!
Thanks. I don’t keep a lot of old Q-Chem versions around, oldest one I have right now is ca. 5.3, and that one assigns zero to the ghost atom radius even with your input file, same as current version Q-Chem 5.4. It looks like this capability (which I wasn’t certain even existed in the first place) was deprecated at some point since 5.0, probably inadvertently. I will submit this as a bug ticket.
By the way - be very careful when using this with TDDFT. Because the ghost atom has basis functions, your “empty” cavity has the ability to support density and given TDDFT’s propensity for low-energy CT excitations, I can definitely see this going awry. I would probably only do this calculation with range-separated hybrids, maybe not even CAM-B3LYP (which doesn’t have quite the right asymptotic behavior for a charge-separate state, because it violates alpha+beta=1 constraint) but rather something like LRC-wPBE or LRC-wPBEh that will properly penalize long-range charge transfer.
I have identified the code changes where this functionality was disabled (to facilitate something else), so it’s a slightly non-trivial fix and I don’t have an estimate, sorry.
Hi, when we put the PCM_PRINT >=2， we don’t get the pqr files, right? the pqr file is in the output file? Could you explain to me how to use VMD to check the cavity? thanks too much.
PQR file is in the Q-Chem output. Starts with a line that reads
----- Tesselation .PQR file -----
and ends with
----- End of Tesselation .PQR file -----
(Just search for PQR.) That PQR file has the grid points, with sizes that are proportional to the weights. Following that, there is a second PQR file that starts with
----- Cavity ESP .PQR file -----
This one contains the electrostatic potential on the cavity surface.
Thanks so much, my current output files can not find this words. the below is key words in my input files.I will try the input files in this question.Thanks so much.
$rem
JOBTYPE EDA
EDA2 1
method wB97X-V
BASIS def2-TZVPD
XC_GRID 000099000590
NL_GRID 1
UNRESTRICTED FALSE
MAX_SCF_CYCLES 200
SYMMETRY FALSE
SYM_IGNORE TRUE
MEM_STATIC 2000
BASIS_LIN_DEP_THRESH 6
THRESH 14
SCF_CONVERGENCE 8
SCF_PRINT_FRGM TRUE
SOLVENT_METHOD PCM
$end
I use Q-Chem 5.2 (devel), Q-Chem, Inc., Pleasanton, CA (2019) to run Ghost_Atoms_3.inp calculations.
the job will stop at the below. Is it the version’s problem?
ATOM 735 GH UNK 11 11 -0.942591 -2.013892 -2.013892 -nan 0.398224
ATOM 736 GH UNK 11 11 2.013892 0.942591 2.013892 -nan 0.398224
----- End of Tesselation .PQR file -----
To get the PQR files, you need to set PCM_PRINT = 2 in the $rem. For the other problem, I cannot diagnose it without your full input file including the $molecule.
