Adding missing parameters for polymers with OPLS-AA using quantum mechanical (QM) calculations

EMC
1

Hi all,

I’m trying to build a polymer structure using EMC with OPLS-AA force field, but I’m encountering missing parameter rules. The build works fine with PCFF, but fails with OPLS-AA. I’ve successfully added custom rules for c3ab before, but I’m not sure about the proper approach for these specific cases.

SMILES structure 
left_part  c1(c4c(*)c[N+]cc4)ccc(c2ccc(c3ccc(*)cc3)cc2)cc1, 1, left_part:1, 2, left_part:2 
cap_H *H, 1, left_part:1, 1, left_part:2

ERROR messages
Warning: increment pair {c3ab, c4} not found.
Warning: increment pair {c4, h1a} not found.
Warning: no rule found for {group, site} = {left_part, 7}.
Warning: no rule found for {group, site} = {left_part, 8}.

Any guidance on the proper syntax for adding these rules or alternative approaches would be greatly appreciated!

Thanks in advance,
Kostas

2

Hi Kostas,

Which version of OPLS-AA do you use? 2012 or 2024? I would assume the latter, judging by your reported error?

3

Hi,

I’m using the latest version of OPLS-AA (2024). Do you have any suggestions for handling these errors?

4

In general, bond increments are used to determine partial charge. Each bond provides the increment of the partial charge of a central atom, thus all adding up to the final partial charge of that central atom.

Specifically, for your case, you would need to know the partial charge of the atom in order to assign the correct bond increment. You could use field_increment warn or ignore to assign a zero value. Though, I would do that only when you are sure this is ok. In your case, I would not expect the c4 h1a contribution to be 0. It should be of the ilk of the c4 h1 contribution. You could check the available bond increments and assign missing ones based on existing ones.

Missing rules will have to be provided for the indicated atoms. You could write out a debug.emc by adding

ITEM EMC stage=0 spot=2
put = {name -> “debug”};
ITEM END

to your .esh. You can then check which atoms exactly are meant by checking out the (* Groups *) section.

5

Hi again,

I have a few follow-up questions since I’m still learning this workflow.

I’m working on poly(aryl piperidinium) structures for anion exchange membrane (AEM) applications. Getting the partial charges correct is critical for my simulations since ionic conductivity and water uptake are highly sensitive to the charge distribution, especially around the quaternary nitrogen cation and the aromatic backbone. That’s why I’m trying to be very careful about the charge assignment procedure.

I have calculated the partial charges for all atoms in my structure (repeating unit) using QM. Is there a way to bypass the bond increment procedure entirely and directly assign these charges?

This might seem basic, but I’m struggling to understand how to map my QM charges to the correct atoms in EMC. Specifically:

  1. How can I identify which atom in my .esh file corresponds to which atom in my QM calculation?
  2. Does EMC provide any output that shows atom indices/IDs?
  3. Is there a way to visualize or export the structure with atom numbers before building, so I can create the correct mapping between my QM calculation and EMC’s atom ordering?
6

You can assign your calculated quantum mechanical (QM) partial charges directly to your SMILES representing your chemistry under ITEM GROUPS. The directive charges controls user-defined partial charges with the two main options additive or override. In the former case, user-defined charges are added to whatever a force field defines, whereas override ignores the force field assigned charges and overrides them with the give user-defined input. The following

ITEM	GROUPS
etoh:charges=override	[C-0.3]([H+0.1])([H+0.1])([H+0.1])[C+0.03]([H+0.1])([H+0.1])[O-0.65][H+0.42]
ITEM	END

serves as an example for redefining the partial charges of ethanol. Note, that you will have to specify your hydrogens explicitly. When omitted, their partial charges will be set to zero in the case override. You can also check ‘Simulation Setup => File Formats => Chemistry File => Groups’ in the EMC manual.

With respect to your listed questions:

  1. I would use a standard SMILES representation in order to associate the correct QM partial charge with the correct atom
  2. By default EMC outputs PDB/PSF, which both include atom IDs
  3. EMC represents atom in memory in the order in which they are provided by your SMILES representation; this order is also reflected in the aforementioned PDB/PSF representation

PDB/PSF can be viewed with VMD. For the user’s convenience, EMC provides a .vmdwhich can be used directly in combination with VMD through its -e option.

7

Hi again,

As I mentioned previously, I am working on poly(aryl piperidinium) structures. Over the last few days, I have been able to create the following structure without problems using OPLS-AA 2024:

right_part c1(c(c4ccccc4)(c(F)(F)F)(*))ccc(c2ccc(c3ccc(*)cc3)cc2)cc1, 1, left_part:2, 1, right_part:2, 2, right_part:1, 2, left_part:1
left_part c1(*)ccc(c2ccc(c3ccc(*)cc3)cc2)cc1, 1, right_part:2, 1, left_part:2, 2, right_part:1, 2, left_part:1
cap_H *H, 1, left_part:1, 1, left_part:2, 1, right_part:1, 1, right_part:2

The problem arises when I try to incorporate a cationic piperidinium group into the left_part using the following format:

left_part c1(C4C(*)C[N+]CC4)ccc(c2ccc(c3ccc(*)cc3)cc2)cc1, 1, right_part:2, 1, left_part:2, 2, right_part:1, 2, left_part:1

The warnings I receive are:

core/field/apply/rules.c:387 match = 0x12596bf10
core/field/apply/rules.c:453 type[45] = h1a
Warning: no rule found for {group, site} = {left_part, 9}.
Warning: no rule found for {group, site} = {left_part, 10}.
Error: core/fields.c:433 FieldsApply:
       Missing rules.
       Program aborted.

I believe the issue is related to the hydrogens bonded to the nitrogen, which is in turn bonded to the carbon in the piperidinium ring. Do you have any suggestions on how to handle this problem? Should I modify the .top or .prm files for OPLS-AA 2024? I’ve successfully added custom rules for c3ab before, but I’m not sure about the proper approach for these specific hydrogen-nitrogen bonding cases.