Maintaining charge balance during system generation with EMC

EMC
1

Hi Everyone,

I am trying to generate random copolymer and terpolymer systems with negatively sulfonate groups and positively charge hydronium groups, however, I am having a difficult time inserting the correct number to hydronium ions to achieve charge balance. Whenever I use EMC to generate random polymer systems, a different number of sulfonate groups are generated. I use the attached .esh file 3 times and there were 52, 48, and 56 sulfonate groups generated each time. This makes it impossible to calculate the correct number of hydronium to insert in my system to maintain charge balance. Is there anyway to generate only polymers first so I can count the number of sulfonate groups and then insert the hydronium ions? I have attached my input script to this post.

Thank you for the help,
polymer_test.esh (1.2 KB)

Max

1 Like
2

Hey Max -

If you include a line at the top of your input file,

field_charge false

then EMC will build your structure without erroring out about the electroneutrality. What you could then do is look at the .data file and count the number of sulfonate groups with a python script, then manually add the hydronium ions.

I would say, build a polymer system with the field_charge false flag and only one hydronium group so its parameters get included in the force field file, then use a Python script to count the number of sulfonate groups S in the system and add S-1 more hydronium ions at random coordinates and rotational angles (and also append the correct topology to the data file). This would allow you to automate the process smoothly.

Cheers,

Sam

3

Hi Max,

You can warrant charge neutrality in one of two ways. The first way is to include the counter ion in your monomer definition, which would change your ITEM GROUPS paragraphs as follows:

# Groups section

ITEM	GROUPS

poly1		*CCC(c1ccc(S(:O)(:O):O)cc1)CC*.O(H)(H)(H), &
			1,poly1:2, 2,poly2:1, 2,poly3:1
poly1_e1	CCC(c1ccc(S(:O)(:O):O)cc1)CC*.O(H)(H)(H), &
			1,poly1:1, 1,poly2:1, 1,poly3:1
poly1_e2	*CCC(c1ccc(S(:O)(:O):O)cc1)CC.O(H)(H)(H), &
			1,poly1:2, 1,poly2:2, 1,poly3:2

poly2		*CCCCC*, &
			2,poly1:1, 1,poly2:2, 2,poly3:1
poly2_e1	CCCCC*, &
			1,poly1:1, 1,poly2:1, 1,poly3:1
poly2_e2	*CCCCC, &
			1,poly1:2, 1,poly2:2, 1,poly3:2

poly3		*CCC(c1ccccc1)CC*, &
			2,poly1:1, 2,poly2:1, 1,poly3:2
poly3_e1	CCC(c1ccccc1)CC*, &
			1,poly1:1, 1,poly2:1, 1,poly3:1
poly3_e2	*CCC(c1ccccc1)CC, &
			1,poly1:2, 1,poly2:2, 1,poly3:2

water		O
hydronium	O(H)(H)(H)

ITEM	END	# GROUPS

The mass of one monomer includes the additional cation, thus allowing for inclusion in the mass fraction calculation. The inclusion method also allows you to set system_charge true under ITEM OPTIONS.

An alternative way is to redefine the number of hydronium clusters after the fact. In this case, you will have to use system_charge false under ITEM OPTIONS. I would define my clusters section as

ITEM	CLUSTERS

poly		random, 10
water		water, 432
hydronium	hydronium, 0

ITEM	END	# CLUSTERS

An extra verbatim EMC section is added to check for the system charge and calculate the correct number of hydronium molecules

ITEM	EMC	phase=1	spot=1

variables	= {
  n_hydronium	-> int(abs(charge())+0.5)
};

clusters	= {
  cluster	-> {
    id		-> hydronium, system -> main, group -> hydronium, n -> n_hydronium}
};
 
field		= {
  mode		-> apply,
  check		-> {
    atomistic	-> true,
    charge	-> true
  },
  debug		-> false
};

ITEM	END	# EMC

The keywords phase and spot control where this verbatim section appears in your resulting build.emc. I added the field typing paragraph with checking system charge neutrality switched on (charge -> true) to check the validity of my calculation. Please note that with this approach the hydronium molecules are excluded when assuming mass fractions for your clusters.

4

Hi veld,

Thank you for your recommendations. I will try them out!

Thanks,
Max