usnistgov / pyPRISM

A framework for conducting polymer reference interaction site model (PRISM) calculations
Other
38 stars 20 forks source link

A confused problem about omega.Intermolecular and omega.NoIntra #22

Closed maojinyuan closed 1 year ago

maojinyuan commented 1 year ago

Hi there,

I'm a new user of pyPRISM, I decided to place two types of polymer(polymer1 N=10 and polymer2 N=150, total density=0.2) and one nanoparticle to a system. Below is my script and the g(r) fucntion. However, I felt something wrong in g(r) function. Here, I just can't understand the difference of omega.Intermolecular and omega.Intra, cause the omega represents the Intra-molecular function.

#!/home/mjy/.miniconda3/bin/python
import pyPRISM
import matplotlib.pyplot as plt
import numpy as np

s_len = 10
l_len = 150
s_n = 5000
l_n = int((0.2*80**3-s_len*s_n)/l_len)

sys = pyPRISM.System(['polymer1','polymer2','particle'],kT=1.0)
sys.domain = pyPRISM.Domain(dr=0.01,length=2048)

sys.density['polymer1']  = np.around(s_n*s_len/(80**3),4)
sys.density['polymer2']  = np.around(l_n*l_len/(80**3),4)
sys.density['particle'] = 1/(80**3)

sys.diameter['polymer1']  = 1.0
sys.diameter['polymer2']  = 1.0
sys.diameter['particle'] = 10.0

sys.omega['polymer1','polymer1'] = pyPRISM.omega.DiscreteKoyama(sigma=1.0,l=1.0,length=s_len,lp=4/3)
sys.omega['polymer1','polymer2'] = pyPRISM.omega.InterMolecular()
#sys.omega['polymer1','particle'] = pyPRISM.omega.InterMolecular()
sys.omega['polymer1','particle'] = pyPRISM.omega.NoIntra()

sys.omega['polymer2','polymer2'] = pyPRISM.omega.DiscreteKoyama(sigma=1.0,l=1.0,length=l_len,lp=4/3)
#sys.omega['polymer2','particle'] = pyPRISM.omega.InterMolecular()
sys.omega['polymer2','particle'] = pyPRISM.omega.NoIntra()

sys.omega['particle','particle'] = pyPRISM.omega.SingleSite()

sys.potential[sys.types,sys.types] = pyPRISM.potential.WeeksChandlerAndersen(epsilon=1.0)
#sys.potential['particle',sys.types] = pyPRISM.potential.WeeksChandlerAndersen(epsilon=1.0)
#sys.potential['particle','particle'] = pyPRISM.potential.WeeksChandlerAndersen(epsilon=1.0)

sys.closure[sys.types,sys.types] = pyPRISM.closure.PercusYevick()
sys.closure['particle','particle'] = pyPRISM.closure.HyperNettedChain()

PRISM = sys.solve()

rdf = pyPRISM.calculate.pair_correlation(PRISM)
plt.plot(sys.domain.r,rdf['particle','polymer1'])

plt.show()

屏幕截图_20221217_143724

tgartner commented 1 year ago

Hi, sorry for the confusion. pyPRISM.omega.InterMolecular() and pyPRISM.omega.NoIntra() are exactly the same, both denoting an intermolecular correlation function of zero everywhere, i.e. sites never appear in the same molecule. The InterMolecular() syntax was from an older version of pyPRISM, and was changed to NoIntra() to avoid exactly the confusion you are having. If you look at https://github.com/usnistgov/pyPRISM/blob/master/pyPRISM/omega/InterMolecular.py, you can see that InterMolecular() is now just an alias of NoIntra(). In terms of the polymer-particle g(r), were you expecting something significantly different? Seems mostly reasonable to me given the relative diameters of polymer and particle sites. I would suggest increasing the length quite a bit when defining sys.domain (perhaps to 8192) so that the domain is large enough for the correlations to decay to one over the length scale of the domain.

maojinyuan commented 1 year ago

Thank you for your response.