Issue with Enhanced Monte Carlo and partial charges

Hello!

I face some issues regarding the automatic application of partial charges with EMC (v9.4.4).

Below, I attach a minimal example of an acid (p. 123 from the manual) modeled with the PCFF force field.

ITEM	OPTIONS
replace			true
number			true
field_charge	false
field			pcff
#field			opls
build_dir		.
project			pani
ITEM	END

ITEM	GROUPS
#acid:c=a		[C](=[O])[O-1] # does not work
acid:c=a		[C](=[O])[O-0.99999999] # works
ITEM	END

ITEM	CLUSTERS
acid			acid,1
ITEM	END

Note that when I set the charge to -1 I receive the error msg:

Info: groups = {ngroups -> 1, group ->
        {id -> acid, chemistry -> "[C](=[O])[O-1]", depth -> 8, charges ->
          additive, charge -> 0, terminator -> false}, ndeletes -> 0,
        npolymers -> 0}
Info: field = {mode -> apply, style -> none, error -> true, debug -> false,
        check -> {atomistic -> true, charge -> false}}
Info: applying '~/Programs/emc/v9.4.4/field/pcff/pcff.frc'
Warning: increment pair {c=2, o-} not found.
Error: core/script/field/entry.c:645 ScriptFieldEntryApply:
       Missing force field parameters.
       Program aborted.

Changing the option from c=a[daptive] or c=o[verwrite] does not address the issue.

On the other hand, setting the charge to -0.999999 works.
The example with the water reported in the manual works as well.

It seems that the program attempts to find a specific rule for {c=, o-} pairs which is not found in the force field; thus it results an error.

I face the same issue with other chemistries as well; e.g., charged [N+1] in emeraldine salts, chlorine with variable charge (e.g., [Cl-0.7], etc.).

The solution of using a charge equal to 0.9999999 does not work for all cases; e.g., [Cl-0.9999999] does not work.

This is not a major issue since one can correct the charges by hand posteriorly, but I’m wandering whether there is a canonical way to do this.

Thanks,
Aris

Aris,

Thank you for your bug report. I think you might be misinterpreting the functionality of the charges option under GROUPS. Let us focus on your acid group definition,

ITEM		GROUPS
acid:c=a		[C](=[O])[O-1]
ITEM		END

For generality’s sake: the letter c can be exchanged with charges; the letter a can be exchanged with additive. This setting will add a charge of -1 to whatever the force field has assigned.

Starting the group assignment from scratch: the following defines deprotonated formic acid can be defined through

ITEM		GROUPS
acid		C(=O)[O-].[Na+]
ITEM		END

Alternatively, using partial bonds through C(:O):O.[Na+] works as well. I normally add the counter ion in the same definition as the acid, thus avoiding later system neutrality issues. I also replaced the charge designation of -1 to - as to avoid internal interpretation of the numerical value. The use of a numerical value can end up in some rounding issues depending on your platform of execution. I would use [Cl-] instead of [Cl-1] as mentioned in your example. Note that EMC transforms C(=O)[O-] to C(:[O-0.5]):[O-0.5] internally.

The use of the charges keyword allows to influence whatever partial charge EMC assigns. Three options are allowed:

  • a or additive: add the provided charges on top of whatever the force field assigns,
  • o or override: use the provided charges instead of the force field assignments,
  • f or field: ignore the provided charges and use the force field assignments instead.

Adapting charges should be handled with care. Normally there is no reason to specify the charges keyword, in which case I would suggest to omit the keyword altogether. Having said that, I created the following for illustrational purposes:

# Options section

ITEM	OPTIONS

replace			true
field			pcff
density			0.1
number			true
system_charge	false
focus			true
pdb_licorice	true
emc_execute		true

ITEM	END	# OPTIONS

# Shorthand section

ITEM	SHORTHAND

acid:charges=field	CC(=O)[O-1].[Na+],1

ITEM	END	# SHORTHAND

The option charges=field represents the default and can be omitted. This generates the following standard charge neutral entry in a .psf (partial charges are listed in column 7),

       8 !NATOM
       1 acid 1    acid C    c             -0.389       12.0112       0
       2 acid 1    acid H1   hc             0.053       1.00797       0
       3 acid 1    acid H2   hc             0.053       1.00797       0
       4 acid 1    acid H3   hc             0.053       1.00797       0
       5 acid 1    acid C1   c-            0.2974       12.0112       0
       6 acid 1    acid O    o-           -0.5337       15.9994       0
       7 acid 1    acid O1   o-           -0.5337       15.9994       0
       8 acid 1    acid Na   na+                1         22.99       0

Changing field to additive results in a .psf with a total charge of -1,

       8 !NATOM
       1 acid 1    acid C    c             -0.389       12.0112       0
       2 acid 1    acid H1   hc             0.053       1.00797       0
       3 acid 1    acid H2   hc             0.053       1.00797       0
       4 acid 1    acid H3   hc             0.053       1.00797       0
       5 acid 1    acid C1   c-            0.2974       12.0112       0
       6 acid 1    acid O    o-           -1.0337       15.9994       0
       7 acid 1    acid O1   o-           -1.0337       15.9994       0
       8 acid 1    acid Na   na+                1         22.99       0

Note the change in charge with -0.5 for both o- types. Using override instead of field generates

       8 !NATOM
       1 acid 1    acid C    c                  0       12.0112       0
       2 acid 1    acid H1   hc                 0       1.00797       0
       3 acid 1    acid H2   hc                 0       1.00797       0
       4 acid 1    acid H3   hc                 0       1.00797       0
       5 acid 1    acid C1   c-                 0       12.0112       0
       6 acid 1    acid O    o-              -0.5       15.9994       0
       7 acid 1    acid O1   o-              -0.5       15.9994       0
       8 acid 1    acid Na   na+                1         22.99       0

Again, note that EMC transforms C(=O)[O-] to C(:[O-0.5]):[O-0.5] internally.

Dear Dr. In 't Veld,

Many thanks for your swift reply. I will certainly take your suggestions into account, and I’m confident they will be effective.

I truly appreciate your continued involvement in the project and your generosity in sharing your insights with the community.

Thanks,
Aris

Dear Aris,

My pleasure. Thank you for your kind words. I hope my exposé helped you understand the difference in functionality.

Best,
Pieter