Questions about CHARMM Drude/SWM4 water treatment in ChemShell POLCOS

dguion · 24 May 2026 00:34

Dear ChemShell developers,

I am setting up a QM/MM calculation in Tcl-ChemShell 3.7.2 using the CHARMM Drude polarizable force field through the documented mm_polcos interface:

mm_polcos=yes
polcos_efield=drude
polcos_atom_polcosq= $cospara

My understanding from the manual is that the POLCOS/Drude model is specified using polcos_atom_polcosq, where each entry gives the polarizable parent atom ID, the polarizability, and the massless Drude/COS charge:

set cospara {
{ atom_index alpha_au q_drude }
…
}

and this is then passed as:

set polarg [ list \
polcos_scale14=1.0 \
polcos_maxcycle=10 \
polcos_toler_energy=1.0e-6 \
polcos_maxdx=2.0e-5 \
polcos_rmsdx=1.0e-5 \
polcos_efield=drude \
polcos_inmaxcyc=1000 \
polcos_atom_polcosq= $cospara \
]

I also inspected the ChemShell 3.7.2 source code to check my interpretation. In src/hybrid/polcos.c, polcos_atom_polcosq appears to be parsed as a three-field list containing the polarized MM atom ID, polarizability, and COS charge:

/* polarized MM atom ID, polarizability and COS charge*/
…
b->atomid = iatm - 1;
b->alpha = alpha;
b->cosq = cosq;

The same file appears to initialize the COS position from the parent atom position:

/* set initial COS position to be polarized MM atom position */
…
a = f->atoms + iatm;
for(j = 0; j < 3; j++) polcosf->atoms[i].pol_pos.x[j] = a->pos.x[j];

For polcos_efield=drude, the code then appears to use the gradient/electric field associated with the COS charge:

/* For Charmm */
…
ix[1] = calc[0].ntot + iatom;
MATRIX_get_value(qmgrad,ix,&edum);
b->efield_qm.x[ix[0]] = -1.0 * edum/b->cosq;

Based on this, my working interpretation is that, in this workflow, the ChemShell/DL_POLY coordinate object contains the real atoms, while the Drude/COS sites are generated and updated internally by POLCOS from the parent-atom positions and parameters in polcos_atom_polcosq. That is, I am not treating the explicit CHARMM Drude particles as ordinary MM atoms in the ChemShell/DL_POLY coordinate object. I am instead relying on polcos_atom_polcosq to define the corresponding polarizable positions.

I am writing because I encountered a specific problem with SWM4-type CHARMM Drude waters.

When I include the SWM4 water oxygens/parents in cospara, ChemShell proceeds into the POLCOS initialization and iterative polarization cycles:

COS model initialization*******

running polcos_init

followed by:

COS_DOUBLE_ITERATION******************

iteration: 1

=MM SCF cycle: 1
MM SCF converges?
Max cos displacement: …
RMS cos displacement: …

However, when waters are included as POLCOS polarizable centers, the resulting water geometries become unphysical. In particular, some water structures stretch or distort badly during the calculation instead of remaining close to sensible SWM4 water geometry. This does not look like a normal optimization displacement and appears specifically tied to treating the SWM4 water sites as POLCOS polarizable centers.

Since polcos.c appears to initialize each COS charge at its parent atom position and then update it from the field instead of reading or preserving the original explicit Drude/auxiliary-site coordinates, I am concerned that the full SWM4 water-specific representation may not be retained when the water oxygen is treated as a generic POLCOS parent site.

As a practical workaround, I have kept the waters in the MM region but excluded SWM4 residues from cospara. In other words:

Waters remain present in the MM PSF/PDB and participate in the MM calculation.
Auxiliary-site charges from the explicit Drude/lone-pair representation are collapsed onto the corresponding real atoms in the ChemShell/DL_POLY real-atom representation.
Waters may still be mobile or frozen depending on my active-atom selection, excluding them from cospara does not by itself freeze them.
SWM4 water oxygens are not included as POLCOS polarizable centers.
Non-water polarizable atoms are still included in cospara.

To prepare the input, I use a small preprocessing script, collapse_drude_for_polcos.py. Its purpose is to convert an explicit CHARMM Drude PSF/PDB into the representation that I understand ChemShell’s POLCOS interface expects. It removes auxiliary Drude/lone-pair sites from the coordinate and topology files passed to ChemShell, remaps atom indices consistently, collapses auxiliary-site charges onto their parent real atoms for the fixed MM charge list, and writes the corresponding ChemShell input lists, including remapped QM and active-atom selections, real-atom coordinates/topology, atom charges/types/residue groups, and cospara entries of the form:

{ parent_atom_index alpha_au abs(q_drude) }

My intention here is to supply ChemShell with real atoms plus polcos_atom_polcosq, since the documented POLCOS interface and the implementation in polcos.c appear to expect polarizable centers to be specified through parent-atom polarizabilities and Drude/COS charges rather than as explicit Drude particles in the coordinate object.

With this treatment, the calculation runs and the POLCOS iterative calculation converges. The SWM4 waters remain as mobile or frozen MM atoms according to my usual active-atom selection, but they are not included as POLCOS polarizable centers. Under this setup I no longer observe the large water distortions seen when the SWM4 oxygen atoms are included in cospara.

My questions are:

Is my interpretation of the intended ChemShell POLCOS/CHARMM Drude workflow correct: real atoms in the ChemShell/DL_POLY coordinate object, with polarizability specified through polcos_atom_polcosq, rather than retaining explicit Drude particles as normal MM atoms?
Does ChemShell’s POLCOS implementation support CHARMM Drude water models such as SWM4?
For SWM4 waters specifically, what is the intended way to represent them in a ChemShell mm_polcos QM/MM calculation? Should they be passed through polcos_atom_polcosq like other polarizable atoms, or do they require a different treatment?
If the current mm_polcos interface does not fully support SWM4-type waters, is treating those waters as fixed-charge MM sites while retaining POLCOS for the non-water polarizable MM atoms a reasonable workaround, or would you recommend a different approximation?
Are there recommended precautions near the QM/MM boundary if water polarization is omitted but protein/substrate polarization is retained?
More generally, are there known limitations in ChemShell 3.7.2 for CHARMM Drude water models under mm_polcos? In particular, could the large/unphysical water distortions I observe when SWM4 water oxygens are included in cospara be related to the POLCOS implementation initializing COS positions from the parent atom positions, rather than reading or preserving the original explicit Drude/auxiliary-site coordinates from the input structure?

I want to make sure I am not misusing the interface or introducing an unintended approximation. The key symptom I observed is that including SWM4 waters as POLCOS centers gives unstable/unphysical water geometries, while excluding them from cospara gives a stable POLCOS calculation.

Thank you for any guidance.