atomistic PRISM calculation, for determination of interaction parameter

[arashelahi]

Hi
I am new user to pyPRISM and integral equation theory. thus my question may look simple, sorry about that in advance.
I have the following questions.
1-I would like to assign an atom to each site in my monomer, rather than coarse-grained simulation, for doing that, what is the appropriate OMEGA and why?
2-In the examples that I saw and the referenced papers, solvent (in my case water) is simulated implicitly, can pyPRISM get explicit Solvent sites, because eventually, I would like to determine interaction parameter between the POLYMER and WATER, so I need at least two components.
3- In atomistic simulation, how interaction parameter between monomer and water can be determined? Is that an average of interaction between all atoms of monomer with water? or what is that?

The reason I would like to do atomistic simulation is, eventually I want to compare my results with experimental results(radial distribution) that is existing for atoms, rather than group of atoms
thank you for your response

[tgartner]

Hi Arash-
Thank you for your questions.

1. PRISM theory can certainly be used at an atomistic resolution, however, you are correct that defining the intra-molecular correlation functions (\omega_i,j) will be a challenge. For your system, there will likely not be a pre-existing \omega_i,j that you will be able to define for each pair of atom types i and j. Instead, you will need to run an atomistic simulation of your system of interest and calculate the \omega_i,j for all pairs of sites from the results of your simulation. Then, you can use these calculated \omega_i,j in pyPRISM by using the pyPRISM.omega.FromFile method. You can refer to section II.C (specifically, equation 12) of our pyPRISM paper for more details about this calculation.
2. Explicit solvent species is also possible, but as with the polymer, you will need to specify the \omega_i,j correctly for all site types to capture the intra-molecular correlations (shape) of the water molecules.
3. This is a good question, as the chi calculation as currently implemented in pyPRISM is strictly valid for a 2-site system. One way to get around this issue is to combine the various direct correlation functions (c_i,j) of the sites within the POLYMER and WATER molecules to calculate direct correlation functions between 'effective' POLYMER and WATER sites. Then, you can use the chi calculation as currently implemented to calculate the chi between these effective POLYMER and WATER sites. Please see the attached pdf file for some theoretical background and more details about this approach. In the future, we plan to add some features to pyPRISM to facilitate this process, as well as add the attached notebook to our pyPRISM tutorial.
   CombiningChi.pdf

[arashelahi]

thank you for your response Thomas, and for the document.
I will work on them and reach you out back.
Just one question.
By "you will need to run an atomistic simulation of your system of interest and calculate the \omega_i,j", did you mean MD simulation?

[tgartner]

Yes, I think a single-chain MD simulation of your polymer of interest solvated in water should be sufficient. You'll just need to make sure that you sample for long enough (get enough independent configurations) to make sure you have good statistics for all pairs of omega_i,j. Good luck!

[arashelahi]

thank you for your response
and I'm finding pyPRISM.trajectory.Debyer module in pyPRISM, can I use the determined trajectories of MD, and let pyPRISM do the calculations?
Also, should I use all the positions in all of the time steps, or some snapshots could be enough. I appreciate if you could explain a bit this part

[martintb]

and I'm finding pyPRISM.trajectory.Debyer module in pyPRISM, can I use the determined trajectories of MD, and let pyPRISM do the calculations?

Absolutely! You should just keep in mind that, while the implementation of the calculation is optimized, it is still extremely computationally intensive. I have another c++ tool (https://github.com/martintb/correlate) that you should consider if you need to apply it to a large molecule.

should I use all the positions in all of the time steps, or some snapshots could be enough.

This is a difficult question to answer. You need to ensure that enough time steps are included to properly sample all of the relevant configurations of the molecule. You should maybe test out including different numbers of snapshots and see how this affects the final omega.

[arashelahi]

thank you Martin
Due to intensity of these calculations, at the same time, I am trying to simulate a coarse-grained polymer water solution, and using pyPRISM.omega.Gaussian and Freelyjointchain modules, but I had convergence problem, can only OMEGA be the reason of divergence?
I am still not sure about the rest of the parameters I set, that are correctly scaled or not. I appreciate if I could share my CG script, and if possible you comment on them.

[tgartner]

Convergence can be affected by many factors, including the spatial discretization (dr) and number of points used in defining the system, the strength of site-site interactions, the concentration of the system, choice of omega, choice of closure, etc. Go ahead to post your code, I can look/comment on it when I get a chance. Please include as much information about your system as possible so we can understand what you are trying to do.

[arashelahi]

thank you Thomas.
I tried to recalculate all the information before posting here, but still I have convergence issue.
I am posting the code here, above each command, the explanation and calculations are provided.
but in general this is our case:
I have a polymer with 40 repeating units of -(C2H4O)- , solvated in water. the sites are monomer and water molecules. explanation about density all provided in thee code. Besides, the data I used as forcefield parameters and distances are all from table 1 of the following paper:
https://pubs.acs.org/doi/abs/10.1021/jp9058966
Although, all the used data provided as comments in the code.
thank you so much
polymer_water.txt

[tgartner]

I have several comments:

1. It looks like you transposed the epsilon and sigma values you input to the LJ potential. In your comments you suggest a water-water sigma of 4.7/3.1, but in your definition of sys.potential['water','water] you use sigma=2.02, which does not equal 4.7/3.1.
2. PRISM may have difficulty converging with such strong attractive interactions, especially starting directly with these high epsilon values. You may wish to start with a weakly attractive system (say all LJ epsilons = 0.1 kT) and then slowly iterate up to higher attraction strengths, using the PRISM output of the previous iteration as the input to the next. For examples on how to structure your code, see the tutorial notebooks NB5 or NB6, which iterate over various system parameters to achieve convergence with the final desired system.
3. Similarly, you may have better success starting with a system of initially higher polymer concentration, then iterating down to a lower concentration.
4. I would double check how you defined the overall density of your system. It looks like you made some calculations based on comparison to experimental values, but your system seems to be at an unphysically high density. Taking the corrected water sigma in the LJ potential, the total occupied volume fraction (JUST accounting for the water) is (volume of one water site)(site density of water), which based on your LJ potential and density settings is ((4/3) pi (1.52/2.0)^3)(1.0) = 1.82. So, there appears to be SIGNIFICANT intermolecular overlap. Usual liquid-like occupied volume fractions with these types of models are more often on the order of 0.4~0.5.

[arashelahi]

Dear Thomas
thank you for your thorough response,
1- that's right. I mistakenly, put epsilon and sigma value reversely, the correct one is sigma=4.7/3.1=1.52, that I put for sigma.
2- thank you so much. I did so, and it got converged (provided following code), but for low density. not high density. I used 0.3 in the following code

3- I will try that, but for now, it could converge with low density of water.
4- I guess I had some misunderstanding about site density, I calculated it in a bad way, density with considering one as reference, but still, using the following way, the the total occupied volume fraction is high.
using d=3.1 Angstrom as ref.
4 chains of polymer, with 40 repeating units and 8000 water molecules in a box of sized=646464 (cubic angstrom)
site density of water=8000/(64/3.1)^3=0.9
site density of polymer=160/(64/3.1)^3=0.018
volume of one water site=((4/3) pi (1.52/2.0)^3)=1.64
Using these values as density, still results in divergence.
However still looks like, with lower density the radial distribution is negative at some points
polymer_water_2.txt

[tgartner]

Good, glad that it is converging in some cases. You can also apply the same iteration idea with the density-- start with the converged state at low density and iterate up to as high a density as you can.

I'm a little confused with the density question overall, however. The bulk density of the water beads makes sense (8000 water molecules in a 64 angstrom cubic box is close-ish to the experimental density of water), but something is strange with the sigma used in the water-water LJ potential. An LJ fluid of number density ~0.9 with an LJ sigma of ~1.5 is far too densely packed for any MD simulation. I would go back to the paper that lists the Martini model parameters and double check that you are interpreting/normalizing all the values correctly.

[arashelahi]

thank you dear Thomas
well, I found those parameters in more than one paper.
table 1 of this paper:
https://pubs.acs.org/doi/pdf/10.1021/jp9058966
first paragraph of section 4.1.2 of this paper: (paper I provided the link already)
https://arxiv.org/pdf/1901.04413.pdf.
martini itp file for this polymer, in the [ nonbond_params ] block of the following link
http://www.cgmartini.nl/images/parameters/polymers/Linear/PEG/martini_v2.0_PEO_PS_CNP.itp.

thank you so much for double checking, I referred to these papers for those values.
Also, would you please briefly explain, why correlation parameters are showing up as negative values, Although it's not diverging. I guess there is no physical meaning for negative radial distribution.

[tgartner]

Ah, I think I got it! Martini water has 4 water molecules per water bead. So, your water site density is 4x higher than it should be. Instead, you should have site density of water=(8000/4)/(64/3.1)^3=0.227.

In terms of the negative g(r), you are correct that this is an unphysical result. There should not be any negative g(r) values. This is a signal that there are numerical issues preventing a good (physical) solution. This can sometimes be solved using the iteration scheme I mentioned above, where you start from a good solution and iterate to your desired system parameters. Or, you can try adjusting the number of points and grid spacing (dr) used in your definition of the system object. This can sometimes affect the numerical result. The key challenge in correctly using PRISM theory is testing these various parameters until you achieve a stable, physical result. Having another method or set of data (either from simulations or experiments) can be helpful in validating your results.

[arashelahi]

Dear Thomas
Thank you for your response
that's right. each water bead includes 4 water molecules in MARTINI force-field. thank you for your following.
Now, for making sure that I am in the right path, I am going to run an MD equivalent system, to compare radial distribution or other properties, resulted from these two methods. after getting some results, I would be back to you.
Another Simple question I have about pyPRISM.
my final goal is to eventually get chi parameter between PEO and Water, and use the resulted chi in Flory-Huggins eqaution for determination of phase-diagram and other characteristics of the system. I know the chi resulted from PRISM theory, is equivalent to SANS interaction parameter, not FLORY parameter. but in pyPRISM documentation, I am seeing it's stated that the package can calculate "FLORY" interaction parameter.
I like to know why we can call it Flory interaction parameter. Are they equivalent in the systems pyPRISM can work on?
thank you so much

[martintb]

So, the difference between the Flory-Huggins chi and the effective chi that pyPRISM calculates is that the PRISM chi includes entropic and enthalpic effects and the Flory-Huggins chi only includes enthalpic effects. I believe we are relatively careful throughout the documentation of always using the phrase 'effective chi' but sometimes we say 'Flory effective chi'.

With that said, the traditional Flory-Huggins chi is simply related to difference in well-depths of your interaction potentials. You should be able to calculate the Flory-Huggins chi from the values you set in the pyPRISM calculation.

Sorry for the delayed response. Cheers! -Tyler

[arashelahi]

Dear Thomas and Martin
thank you for your responses and helps so far, and explanation for chi in the above comment was so clear. thanks
In the past couple of weeks, I had to develop a good-working MD model for PEO-WATER and another polymers solutions, that I can compare them with my PRISM calculation too. now for PEO-WATER, my model is working pretty good, and can validate most of the results of the following paper (including figure 2 and table 3 results).
https://pubs.acs.org/doi/pdf/10.1021/jp9058966
the picture attached below is the comparison between rdf calculated from MD and PRISM theory, for a system of PEO with 36 mers in a 10 nm per side box filled with water. the trend of g(r) is in an agreement, but still there is a considerable deviation in the graphs. I have the following questions.

1- I like to know your opinion about this result. Is that what we can call agreement of MD and PRISM?
2-does using trajectory of MD simulation work for this goal, instead of using FreelyJointedchain form factor?

3- is that possible to provide g(r) too in pyPRISM determined from MD as an input, and instead of iterative calculation of h(r) and c(r), only calculate c(r). I guess this way I could be closer to my MD simulation environment. do you think that could work?

in the following my code, including information about the problem is provided. although, this could not converge by only one kicking off, I needed to cut the energy iteration counters to more intervals, but this is the my main script.

P_W_36mers.txt
.
thank you so much for your following up. I really appreciate it

[tgartner]

You may be able to improve the agreement between MD and PRISM if you try other closures, though I have had reasonable success with the PY closure with similar systems. Additionally, since you have already run MD simulations with these systems, you could calculate the chain omega from the MD results and use the calculated omega within your pyPRISM calculation using the Omega.FromFile method, rather than using the ideal freely jointed chain model (see for example this paper: https://pubs.acs.org/doi/10.1021/acs.macromol.9b00600). The use of this more realistic omega may improve the agreement between the methods. In terms of your other questions, I’m not sure if I understand what you mean by ‘this goal’. Is your goal still to calculate a polymer-water chi parameter? If so, you are correct that you can use MD results to directly calculate h(r) and then solve the PRISM equation to get c(r) and the corresponding chi parameter. However, you don’t need the pyPRISM software for this calculation, as you can code it yourself in a fairly simple script (no iterative numerical solution is required with this approach). See this paper for an example: https://pubs.rsc.org/en/content/articlelanding/2018/sm/c7sm02199b/unauth#!divAbstract However, one downside of this MD-only approach is that very large simulation boxes are needed to provide reliable sampling of long length scale structural correlations that enable the chi calculation. I have found that for dilute polymer solutions, trying to get the c(r) directly from MD was extremely computationally demanding and resulted in large statistical uncertainty. Hope this helps. Best, Tom

…

On Jun 14, 2019, at 1:09 AM, arash303 ***@***.***> wrote: Dear Thomas and Martin thank you for your responses and helps so far, and explanation for chi in the above comment was so clear. thanks In the past couple of weeks, I had to develop a good-working MD model for PEO-WATER and another polymers solutions, that I can compare them with my PRISM calculation too. now for PEO-WATER, my model is working pretty good, and can validate most of the results of the following paper (including figure 2 and table 3 results). https://pubs.acs.org/doi/pdf/10.1021/jp9058966 the picture attached below is the comparison between rdf calculated from MD and PRISM theory. the trend of g(r) is in an agreement, but still there is a considerable deviation. how can I make my PRISM model work better. now I have the following questions 1-does using trajectory of MD simulation work for this goal? 2- is that possible to provide g(r) too in pyPRISM determined from MD as an input, and instead of iterative calculation for h(r) and c(r), only calculate c(r). I guess this way I could be closer to my MD simulation environment. do you think is that a good idea? in the following my code, including information about the problem is provided. although, this could not converge but only one kicking off, I needed to cut the energy iteration counters to more intervals, but this is the my main script P_W_36mers.txt . thank you so much for your following up — You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub, or mute the thread.

[arashelahi]

Dear Thomas
that is correct, "my goal" is still to calculate chi.
However, in this step, I like to integrate MD and PRISM, to figure out if I can replicate existing empirical chi for certain systems, if so, I would shift to pure PRISM calculation, because the the other objective is avoiding extensive MD calculations.
Thanks for your recommendations.
I will try different closures, and also integration of MD and PRISM to compare the results.
thank you so much