I am currently attempting to use EMC to build up a polyimide chain with PCFF force field using partial charges on the atoms that I have derived from quantum mechanical calculations (HF/6-31G* optimization in Gaussian and RESP fitting with Chargemol). To use these charges in EMC, I would like to use the “charge=override” or “c=o” keyword in my groups section where I define the chemistry with SMILES strings, as was shown here. As you can imagine, this is challenging for a monomer and end groups which in my case have ~67 atoms; however, I have done it with some great effort (in the future, I would like to develop a tool to automate this process).
However, I have discovered a bug. despite putting all the atoms, including the now explicit hydrogens, in brackets with their charges (i.e. [C-0.1] or [H+0.1]), the charge for the atom in the SMILES string which immediately follows the polymer connecting point * is NOT read by EMC. Instead, its charge is applied to the following atom, et cetera, until the end of the molecule, and the charge that is supposed to be applied to the final atom is never used. I discovered this using the
field_charge false
flag in my input file and cross-referencing molecule EMC gave me (a monomer of the two end groups merged) versus the optimized geometry and charges from Chargemol.
When I try to assign a charge to the bonding site (i.e., by using ([*+0.0]) in my SMILES string), I get the following error:
Source connect $end1 out of bounds in group
I wanted to flag this because for anyone who wants to do polymeric systems with externally derived partial charges this is a big problem.
Pieter, do you know what the first steps to addressing this bug would be? I would be happy to try and help with it. I’m attaching the relevant files including a slightly modified version of PCFF that just adds an additional rule to type the imide nitrogen as type nb. The .xyz file is my QM optimized structure with RESP charges and the .data file is the EMC output which is not charge neutral due to the bug. DDEC6_even_tempered_net_atomic_charges.xyz (49.7 KB) polymer.data (59.7 KB) polymer.esh (3.1 KB) pcff.frc (332.4 KB) pcff_templates.dat (39.1 KB)
Thank you for reporting this bug. It turned out that – while populating group site charges – I forgot to neglect the asterisks. I updated the code accordingly and successfully checked its validity with the following minimal problem:
# Options section
ITEM OPTIONS
replace true
field charmm/c36a/cgenff
number true
density 0.1
focus true
emc_execute true
ITEM END # OPTIONS
# Groups section
ITEM GROUPS
mono:c=o *[C-0.2]([H+0.1])([H+0.1])[C-0.2]([H+0.1])([H+0.1])*, &
1,mono:2
term:c=o *[C-0.3]([H+0.1])([H+0.1])([H+0.1]), &
1,mono:1, 1,mono:2
ITEM END # GROUPS
# Clusters section
ITEM CLUSTERS
polymer alternate,1
ITEM END # CLUSTERS
# Polymers section
ITEM POLYMERS
polymer
1 mono,2, term,2
ITEM END # POLYMERS
Although the issue was resolved in the previous EMC version (20240801), it appears that a related bug exists in the latest release (20250801). I attempted to generate a benzene ring with a custom charge distribution along with two methyl trimers using the pcff force field:
#=======================================================
ITEM OPTIONS
replace true
field pcff
number true
field_charge false
density 0.05
focus true
emc_execute true
ITEM END
#=======================================================
# Groups
ITEM GROUPS
TP:c=o *[c+0.1]1[c+0.15] ([H+0.05])[c+0.15] ([H+0.05])[c+0.1]*([c+0.15] ([H+0.05])[c+0.15]1([H+0.05])) ,1,TP:2
methyl:c=o *[C-0.8] ([H+0.05])([H+0.05])([H+0.05]),1,TP:1,1,TP:2
ITEM END
#=======================================================
# Clusters
ITEM CLUSTERS
poly random, 1
ITEM END
#=======================================================
# Polymers
ITEM POLYMERS
poly
1 TP,1,methyl,2
ITEM END
#=======================================================
The code is designed to assign artificially positive charges to all atoms. However, in the resulting .psf output file, three atoms (one carbon and two hydrogens) are assigned zero charge. I tested various configurations and found that this issue consistently occurs only in ring structures.
PSF NAMD CMAP
2 !NTITLE
REMARK created by EMC v9.4.4, build Aug 2 2025 16:30:53
REMARK created on Wed Aug 6 17:39:03 2025
18 !NATOM
1 poly 1 TP C cp 0.1 12.0112 0
2 poly 1 TP C1 cp 0.15 12.0112 0
3 poly 1 TP H11 hc 0.05 1.00797 0
4 poly 1 TP C2 cp 0.15 12.0112 0
5 poly 1 TP H21 hc 0.05 1.00797 0
6 poly 1 TP C3 cp 0.1 12.0112 0
7 poly 1 TP C4 cp 0 12.0112 0
8 poly 1 TP H41 hc 0 1.00797 0
9 poly 1 TP C5 cp 0.15 12.0112 0
10 poly 1 TP H51 hc 0 1.00797 0
11 poly 2 meth C c -0.8 12.0112 0
12 poly 2 meth H1 hc 0.05 1.00797 0
13 poly 2 meth H2 hc 0.05 1.00797 0
14 poly 2 meth H3 hc 0.05 1.00797 0
15 poly 3 meth C c -0.8 12.0112 0
16 poly 3 meth H1 hc 0.05 1.00797 0
17 poly 3 meth H2 hc 0.05 1.00797 0
18 poly 3 meth H3 hc 0.05 1.00797 0
The current version of EMC works correctly. You actually made a slight mistake in your TP SMILES. You forgot to put parenthesis around the second connection point:
Leaving out parenthesis means that the connection point is part of the backbone, instead of being a side chain – as it should be. Normally, connection points should not be part of the molecule backbone. With this correction, your full setup.esh becomes
#!/usr/bin/env emc.pl
# Options
ITEM OPTIONS
replace true
field pcff
number true
field_charge false
density 0.1
focus true
emc_execute true
ITEM END
# Groups
ITEM GROUPS
TP:c=o *[c+0.1]1[c+0.15]([H+0.05])[c+0.15]([H+0.05])[c+0.1](*)([c+0.15]([H+0.05])[c+0.15]1([H+0.05])),1,TP:2
methyl:c=o *[C-0.8]([H+0.05])([H+0.05])([H+0.05]),1,TP:1,1,TP:2
ITEM END
# Clusters
ITEM CLUSTERS
poly random, 1
ITEM END
# Polymers
ITEM POLYMERS
poly
1 TP,1,methyl,2
ITEM END
Partial charge assignment occurs as intended with the above, as is illustrated by the following setup.psf snippet:
PSF NAMD CMAP
2 !NTITLE
REMARK created by EMC v9.4.4, build Aug 5 2025 12:56:09
REMARK created on Thu Aug 7 11:07:28 2025
18 !NATOM
1 poly 1 TP C cp 0.1 12.0112 0
2 poly 1 TP C1 cp 0.15 12.0112 0
3 poly 1 TP H11 hc 0.05 1.00797 0
4 poly 1 TP C2 cp 0.15 12.0112 0
5 poly 1 TP H21 hc 0.05 1.00797 0
6 poly 1 TP C3 cp 0.1 12.0112 0
7 poly 1 TP C4 cp 0.15 12.0112 0
8 poly 1 TP H41 hc 0.05 1.00797 0
9 poly 1 TP C5 cp 0.15 12.0112 0
10 poly 1 TP H51 hc 0.05 1.00797 0
11 poly 2 meth C c -0.8 12.0112 0
12 poly 2 meth H1 hc 0.05 1.00797 0
13 poly 2 meth H2 hc 0.05 1.00797 0
14 poly 2 meth H3 hc 0.05 1.00797 0
15 poly 3 meth C c -0.8 12.0112 0
16 poly 3 meth H1 hc 0.05 1.00797 0
17 poly 3 meth H2 hc 0.05 1.00797 0
18 poly 3 meth H3 hc 0.05 1.00797 0
I will add a check to EMC to disallow backbone connection points.