Skip navigation
5800 Views 25 Replies Latest reply: Jan 8, 2014 8:31 AM by reiniera RSS
Currently Being Moderated

Feb 17, 2012 9:08 AM

Script to create Mesocite forcefields using Boltzmann Inversion

Hi everyone,

 

Attached is a script to create Mesocite forcefields using the Boltzmann Inversion method. This method is based on fitting the coarse-grained interactions to spatial distributions, as measured in a fine-grained (i.e. molecular dynamics) trajectory.

 

The input to the script is a coarse-grained 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 coarse-grained to a ring of 4 beads. The script generates a forcefield for coarse-grained 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

Attachments:
  • debelan1 Active Contributor 17 posts since
    Jun 29, 2010

    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

    • jdejoannis Accelrys 565 posts since
      Aug 13, 2007

      Coarse grained forcefields often don't have torsions

      • debelan1 Active Contributor 17 posts since
        Jun 29, 2010

        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

      • debelan1 Active Contributor 17 posts since
        Jun 29, 2010

        The reference that I previously mentioned is "A Coarse-Grained Model for PEO and PEG: Conformation and Hydrodynamics", Lee, de Vries, Marrink, and Pastor, JPCB 2009, 113, 13194.

        Anthony

  • sebastian.tesarski Contributor 5 posts since
    Feb 26, 2009

    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. Intra-intermolecular 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

      • sebastian.tesarski Contributor 5 posts since
        Feb 26, 2009

        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

  • gttmtd21121 Contributor 6 posts since
    Jan 3, 2013

    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

    • jdejoannis Accelrys 565 posts since
      Aug 13, 2007

      Hi Martin,

       

      I have seen this sometimes as well. In my case, there was a lot of scatter in the angle probability distribution P:

      P.PNG

      I also got a PMF (target) curve that does not look right (normally it would have a well around the most common angle):

      PMF (target).PNG

       

      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):

       

      PMF fit.bmp

      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

  • HY Active Contributor 16 posts since
    Apr 26, 2012

    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

      • HY Active Contributor 16 posts since
        Apr 26, 2012

        Reinier,

         

        Thanks a lot! After assigning all the beads to stry, it does work. I also tried the Coarse-graining 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?

      • HY Active Contributor 16 posts since
        Apr 26, 2012

        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 : compute-0-71.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   

        • HY Active Contributor 16 posts since
          Apr 26, 2012

          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 medium-sized 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

           

           

           

            • HY Active Contributor 16 posts since
              Apr 26, 2012

              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 " re-iterate"  ? How should I do the reiteration? Thanks a lot!.

               

              Leon  

  • jasonlee Contributor 3 posts since
    Jan 5, 2014

    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?bond.jpg

      • jdejoannis Accelrys 565 posts since
        Aug 13, 2007

        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.

      • jasonlee Contributor 3 posts since
        Jan 5, 2014

        Hi Reinier,

         

        Thanks a lot! Do you mean the distribution that the result showed was the target All-atom distribution,not the Coarse-grained 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

        • jdejoannis Accelrys 565 posts since
          Aug 13, 2007

          Hi Jason,

           

          Yes the distribution is coming from the all-atom 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.

More Like This

  • Retrieving data ...

Bookmarked By (1)