Feb 17, 2012 9:08 AM
Script to create Mesocite forcefields using Boltzmann Inversion

Like (1)
Hi everyone,
Attached is a script to create Mesocite forcefields using the Boltzmann Inversion method. This method is based on fitting the coarsegrained interactions to spatial distributions, as measured in a finegrained (i.e. molecular dynamics) trajectory.
The input to the script is a coarsegrained trajectory, which can be obtained from a ordinary molecular dynamics trajectory using the TrajectoryCoarseGrainer script, postedhere. The output of the script is a forcefield, listing all the bond, angle, and van der Waals interactions relevant to the input structure. The forcefield parameters are fitted, such that there is a reasonable match with the spatial distributions in the input trajectory. All is documented in a PDF file.
An example is provided for OMCTS, a cyclic tetramer of dimethylsiloxane, which is coarsegrained to a ring of 4 beads. The script generates a forcefield for coarsegrained OMCTS, which can be used in Mesocite for further refinement.
Please note the disclaimer: This custom script is compatible with Materials Studio version 6.0. It is provided "as is" and is NOT supported by Accelrys nor is it warranted for any purpose whatsoever. The user assumes responsibility for any malfunctions or bugs.
Hope it is of use to anyone!
Reinier
Hi Reinier,
Thank you for this very useful and well written script. I noticed that you have incorporated a placeholder of sorts for a subsequent enhancement to include torsional interactions. Have you had the opportunity to implement this? If not, is it something that you are commited to doing?
Regards,
Anthony
Coarse grained forcefields often don't have torsions
Hi Anthony,
Thanks for your interest in this script. You are quite right, this functionality is missing, and still not implemented. Whilst torsion interaction are often ignored in coarse grained force fields, it would still be of interest to validate this for a system of interest.
One line of validation is to measure the torsion distribution in both atomistic and mesoscale, even with a forcefield containing just bond, angle, and van der Waals interactions. When they differ significantly, and reiterating the inversion does not help, you may have to introduce torsion interactions as well. It would be good to have this evidence before extending the script.
You could extend the script along the same lines as for Distances and Angles, by introducing, and implementing, the functions
PrepareTorsion GetFitDataTorsion EstimateParametersTorsion
The first should create torsions monitors in the structure for all relevant combinations of bead types. The second sets up the functional form of interest. The last one requires some thought, as its job is to come up with a reasonable starting point for the Perl routine Algorithm::CurveFit
. From experience, this routine does not always converge without a reasonable starting point.
Best,
Reinier
Hi Reinier (and Jason),
I realize that torsion interactions are typically neglected in coarse grained models but my interest in using them comes from a publication by Pastor et al, in which a torsional interaction was introduced to reproduce the corresponding distributions (and structural characteristics) for PEO and PEG. I don't have the exact reference handy at the moment, but I will post it early next week. So, at least in select cases, there seems to be some evidence that introducing these interactions can be beneficial.
Best regards,
Anthony
The reference that I previously mentioned is "A CoarseGrained Model for PEO and PEG: Conformation and Hydrodynamics", Lee, de Vries, Marrink, and Pastor, JPCB 2009, 113, 13194.
Anthony
Hi Reinier,
I used your script to create a forcefield, but during the process its stops and this is written: "Radial distribution function analysis failed. Intraintermolecular RDF components cannot be calculated for an infinite fragment in Mesocite.Analysis (function/property "RadialDistributionFunction") at e line 611." I tried to change the values of "RDFBinWidth" and "RDFCutoff" but with no luck.
Any ideas where I should look to find a solution?
Kind regards,
Sebastian
Hi Sebastian,
The script is set up to fit the intermolecular radial distribution function, but Mesocite is unable to calculate this for an infinite fragment (e.g. a crystal, or polymer network); it can only do the total (= inter+intra) RDF in that case.
You could experiment fitting the total RDF instead by changing the script
RDFComputeMolecularComponents => "No"; # was "Yes"
RDFChartAsStudyTable>Sheets("Total"); # was "Inter"
No idea if this works as I only tried molecular systems.
Reinier
Hi Reinier,
It works  as for the scripts, but I'm not sure are results are valid.....
You wrote that scripts make only one iteration and probably for better accuracy more iterations are needed. In which part loop should be added? In sub Run? or in other sub?
Sebastian
Hi Sebastian,
One way to iterate this is to run Mesocite with the force field created by the script, and measure the distribution of bonds/angles/nonbonds again. By comparing these distributions with the target ones you can correct the force field, and repeat the process, until convergence.
For example, let u_n(r) be the van der Waals potential of iteration n at distance r, and g_n(r) the corresponding radial distribution function, obtained from a subsequent Mesocite run. Let g_exp(r) be the target distribution. If g_n(r)>g_exp(r) the number of pairs with distance r is too high, and the potential must be increased to correct for this, and vice versa. Hence
u_{n+1} = u_n + C ln(g_n/g_exp)
where C > 0 is a constant, is expected to result in a distribution g_{n+1} that is closer to g_exp.
You can follow this procedure for bonds and angles in the same way.
As only parametrized functions are allowed in Mesocite, you will need to fit u_{n+1}, before you can run the next iteration. This can be done using the FitParameters routine. You can use C to tune the convergence rate.
As a word of warning, this refitting could impact convergence. In practise I found the Perl package Algorithm::CurveFit to diverge occasionally when trying out the above procedure.
Good luck,
Reinier
Hi Reinier,
I’m Martin, I have some question about Boltzmann Inversion code, would you please help me solve the error when I used. The error was happened during the process and it was shown as: "Parameter T0 (193.752) out of bounds, valid range is [0,180] (writing property "T0") at e line 895."
I think that the value of T0 statistics should be at 0~180 degrees, why the T0 was more than 180 degrees when the results were calculated with iteration method via weighting? How can I change the range of parameter T0 and correct the T0 in a range between 0 to 180 degree?
Thanks
Hi Martin,
I have seen this sometimes as well. In my case, there was a lot of scatter in the angle probability distribution P:
I also got a PMF (target) curve that does not look right (normally it would have a well around the most common angle):
The solution was to increase the angle binning from 0.5 to 1.0 so there is less scatter in the angle distribution. Doing that I get a much more reasonable result. Comparing PMF (target) to PMF (fit):
It yields a minimum at 165 degrees.
The error message you got probably comes when trying to put a value that is out of the allowed range (0,180) into the forcefield OFF document.
Best,
/Jason
Hi Reinier,
Thanks a lot ofr writing the script, which is exactly what I need.
I tried your script on a polystyrene polymer structure. The TrajectoryCoarseGrainer.pl script works fine to generate a coarse grained structure with three bead type styrene, styrene(2) and styrene(3). Than I ran the BoltzmannInversion.pl sscript, but the output .off file only contains one type called "styr". Using this .off force fild file, I was unable to run mesocite calculations on the already coarse grained structure with an error msg of "The model contains particles that have not been assigned a forcefield type", which is expected because the coarse grained structure contains three bead types.
Do you know why the script does not recognize all the three bead types? Is this because the names of the three bead types all start with "styr" and somehow the script truncates only the first four letters? Look forward to learning how to fix this. Thanks so much!
Best,
Leon
Hi Leon,
The error message suggests that no forcefield types are set on the bead model. If you select "styr" for all beads you may find it works (Modules  Mesocite  Calculation  Energy  More...  Assign).
Although a bit accidental, it actually makes sense to assign the same forcefield type to different bead types in this case, since they differ only by a hydrogen atom.
However, if you do want to distinguish the bead types in the forcefield, you could set forcefield types on the bead model manually (using the Property Explorer), before running the inversion script. Otherwise, the script will derive them from the bead type names, truncated to 4 characters (a limitation of the forcefield document), which in this case results in the same name.
Alternatively you could use the script in the zip attached to Coarsegraining a trajectory, which has a few extra lines replacing the bead type names with A,B,C,...
Good luck,
Reinier
Reinier,
Thanks a lot! After assigning all the beads to stry, it does work. I also tried the Coarsegraining script you mentioned, it also works great! Now I have another question. I probably can test by myself, but if you already know it will be great. I want to calculate diffusivity of a gas in the polymer. After I coarse grained the structure, the free volume is changed, right? Will that affect the accuracy of diffusivity?
Hi, Reinier,
I ran a mesocite simulation using the obtained force field with a single bead type styr. But the density it predicted is totally wrong, it ends up with 0.027 g/cm3.
I also ran a mesocite simulation (NPT, 300K, 1 bar, 10 ns) with the example files "bulk OMCTS CG " in the folder using its corresponding force field. I found the predicted density is also not correct. It also approaches to 0.059 g/cm3 density. Am I missing something here? Is the obatained force field supposed to work as the atomic force field COMPASS to be able to predict density and other properties? I think I might be missing something. Could you share your thoughts? Thanks a lot.
The output of the bulk OMCTS CG run is as follows:
Mesocite

Task : Dynamics
Version : 6.1
Build date : Oct 12 2012
Host : compute071.local
Threads : Parallel 12
Operating system : linux
Task started : Wed Apr 3 11:53:04 2013
 Dynamics parameters 
Ensemble : NPT
Temperature : 298.00 K
Control method : Nose
Q ratio : 1.0000000
Pressure : 0.00010 GPa
Control method : Berendsen
Decay constant : 1.0000000 ps
Timestep : 1.00 fs
Number of steps : 10000000
Duration : 10000 ps
Initial velocities : Random
Random number seed : 1365007986
 Energy parameters 
Forcefield : \bulk OMCTS CG
Electrostatic terms:
Summation method : Ewald
Accuracy : 0.001 kcal/mol
Buffer width : 0.5 A
van der Waals terms:
Summation method : Bead based
Truncation method : Cubic spline
Cutoff distance : 12.5 A
Spline width : 1 A
Buffer width : 0.5 A
ã€€
 Thermostat 
Nose mass parameter, Q mass: 708.848 kcal/mol*ps^2
ã€€
ã€€
 Dynamics summary 
Initial Final Average Std. Dev.
   
Tot. energy (kcal/mol) 274.478 480.014 438.130 51.253
Pot. energy (kcal/mol) 79.946 121.986 83.708 49.100
Kin. energy (kcal/mol) 354.424 358.028 354.422 14.531
Tot. enthalpy (kcal/mol) 1237.895 506.066 476.823 103.830
Temperature (K) 298.000 301.030 297.998 12.218
Pressure (GPa) 0.130 0.000 0.001 0.003
Volume (A^3) 51521.145 841267.092 409544.477 226328.379
Density (g/cm^3) 0.956 0.059 0.185 0.140
ã€€
ã€€
ã€€
Task terminated : Wed Apr 3 12:34:19 2013
Total CPU time used by Mesocite : 41:08 minutes (2467.83s)
Termination status : Normal
Leon
Hi, Reinier,
To further test the script, I created a structure containing 500 benzene molecules with a density of 0.9 g/cm3 using the amorphous cell tool in MS. Then I ran NPT dynamics for 10 ps with time step of 1 fs, and saved totally 250 frames. Then I corase grained the structures by treating benzene molecule as one bead. The system has only van der waals interaction between the benzene bead. Using the BoltzmannInversion script, I got D0 = 0.160433 kcal/mol and R0 = 5.603 A in the obtained force field file. The D0 value is different from 0.9634 kcal/mol which I found in a paper (Title: "How to predict diffusion of mediumsized molecules in polymer matrices. From atomistic to coarse grain simulations", the authors are Alfonso Gautieri, Simone Vesentini, Alberto Redaelli. If you google it, you can find the full text). I tested the force field and believe D0=0.9634 kcal/mol is the right value. R0 in their paper is 5.62 A, very similar as the one from the script.
Reinier, could you please let me know what could be wrong? If it is possible, could you spend some time to work out the benzene system as well and see what values of D0 and R0 you got? Thanks a lot!!
Leon
I suspect this is a limitation of the Boltzmann Inversion method, and not something you are doing wrong. Basically by fitting the structure you are not guaranteed to get the pressure right (see e.g. http://jcp.aip.org/resource/1/jcpsa6/v114/i2/p1020_s1). Some people make a "pressure correction", but I don't know about the details; you will need to investigate yourself. What is important is to compare the predicted with the target structure (rdfs etc), as you may have to reiterate to get selfconsistency, which is not implemented by the script.
Reinier,
Thanks for your opinion. I had thought the force field obtained should be much better. Is there any requirement for the atomic traj? For example, the total simulation time and intervals between two frames? I suspect that could also affect the resulted force field. Can you talk a litle more regarding " reiterate" ? How should I do the reiteration? Thanks a lot!.
Leon
Hi Leon,
Perhaps the force field is best used in the same ensemble as in which it was obtained. I.e. if you derive NVTpotentials, you should not use them in NPT, without scaling the pressure (running NVT will provide the effective pressure).
Regarding iterations, a simple scheme is
u_{i+1} = u_i − kT ln (g_target/g_i)
where i = 0, 1, 2, ..., g_target is the target distribution and g_i is the distribution as obtained for potential u_i. Starting with an ideal gas (u_0 = 0, g_0 = 1),
u_1 = − kT ln g_target
u_2 = u_1 − kT ln g_target/g_1
u_3 = u_2 − kT ln g_target/g_2
etc. So at distances r where g(r) is too high compared to the target, the potential is increased, and vice versa. As an increase in potential at r will lower g(r), this is expected to converge.
The current script only provides (an approximation to) u_1. You could extend it to do more iterations, but be aware that the data must be fitted to a fixed functional form before you can use it in a next simulation (the script already does this for u_1). Depending on the functional form this may or may not converge.
Best,
Reinier
Hi, Reinier,
I created a structure containing 10 PMMA molecules using the amorphous cell tool in MS. Then I ran NPT dynamics for 500 ps with time step of 1 fs. One PMMA monomer was treated as one bead. Then run your script. The result showed the bond distribtion. But the bond style was harmonic, the distribution must be wrong, How to change the functional form?
Hi Jason,
The graph shows the distribution of the distance between two monomers in the input trajectory, i.e. the one derived from molecular dynamics you are trying to match. It is independent of the functional form you are using to fit it.
However, clearly a bimodal distribution like that will be hard to fit with a harmonic potential. You could try making the beads smaller and see if this gives a distribution with a single maximum.
The script only supports fitting to harmonic functions at present. In principle you could extend it to use other forms, e.g. a quartic potential. I have not tried this, but I suspect you will need to come up with good initial parameter values (in EstimateParametersBondStretch) in order for Algorithm::CurveFit to converge.
Best,
Reinier
Seeing as these are simple PMMA chains with one bead per monomer, I tend to doubt a bimodal distribution is right for the bead connector lengths. Perhaps you want to double check that.
Hi Reinier,
Thanks a lot! Do you mean the distribution that the result showed was the target Allatom distribution,not the Coarsegrained fit distribution? Can I treat all monomers of PMMA as one kind of beads by deleting a hydrogen atom of the tail monomer when building the PMMA model? How to draw the target and fit distribution in a chart just like figure 3 in your documentation using MS?
Many thanks.
Best,
Jason
Hi Jason,
Yes the distribution is coming from the allatom simulation essentially. You can easily reproduce it using the Mesocite  Analysis dialog and the Length distribution task, using the Connector lengths as the objects to analyze. As you can see, it will analyze all connectors regardless of type.
Hi Jason,
I suspect the two maxima reflect the torsion states in the PMMA backbone, e.g. transgauche, transtrans.
To treat all monomers alike you could assign a common force field type to the beads prior to running the inversion script. You can do so via the Property Explorer (Filter = Bead, ForcefieldType = A).
If no force field types are assigned, the script will use the bead type, which is different for head/tail and interior monomers (due to the extra hydrogen atom), and you end up with unnecessary force field interactions.
To create figure 3, I ran Mesocite with the force field resulting from the script. I used Mesocite Analysis to determine the distributions, and plotted them together with those derived from MD. In principle, you could use the difference between the two to correct the force field, and reiterate until convergence. This, however, is not supported in the current version of the script.
Reinier