API tools faq
paste
Advertisement
CasualGamer

COPHY_Main

CasualGamer
Jan 10th, 2021
1,250
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 17.25 KB | None | 0 0
raw download clone embed print report
  1. from __future__ import print_function
  2. import numpy as np
  3. import datetime
  4. import espressomd
  5. from espressomd.pair_criteria import DistanceCriterion #required for cluster analysis
  6. from espressomd.observables import ParticleVelocities
  7. from espressomd.cluster_analysis import ClusterStructure #for cluster analysis
  8. from espressomd.accumulators import Correlator
  9. import espressomd.magnetostatics as magnetostatics
  10. espressomd.assert_features('DIPOLES', 'LENNARD_JONES', 'ROTATION')
  11.  
  12. import tidynamics
  13. import matplotlib.pyplot as plt
  14.  
  15. #helper functions
  16. from functions import combinationRuleEpsilon
  17. from functions import combinationRuleSigma
  18. from functions import averageClusterSize
  19. from functions import addParticles
  20.  
  21. from vtf_files import _writevsf
  22. from vtf_files import _writevcf
  23.  
  24. from scipy.optimize import curve_fit
  25.  
  26. def func(x, a, b, c, d, e):
  27.     return a * np.exp(-b * x)  + c * np.sin(d * x) + e
  28.  
  29. class ParticleType:
  30.     def __init__(self, type, n, density, dipoleMoment, sigma, epsilon, cut):
  31.         self.type = type
  32.         self.n = n
  33.         self.density = density
  34.         self.dipoleMoment = dipoleMoment
  35.         self.sigma = sigma
  36.         self.epsilon = epsilon
  37.         self.cut = cut
  38.     def __str__(self):
  39.         return "type: {}\tn: {}\tdensity: {}\tdipoleMoment: {}\tsigma: {}\tepsilon: {}\tcut: {:.4f}" \
  40.                .format(self.type,self.n,self.density,self.dipoleMoment,self.sigma,self.epsilon,self.cut)
  41.  
  42. def main():
  43.     ################### Parameters start ###################
  44.     # Simulation parameters
  45.     densitySmall = 0.001
  46.     densityLarge = 0.02
  47.     sigmaSmall = 1.  #Diameter of small particle | Definition needed here since multiple ParticleType members depend on it | * 10^-9m
  48.     sigmaLarge = 2.  #Diameter of large particle | Definition needed here since multiple ParticleType members depend on it | * 10^-9m
  49.     dipoleMomentSmall = 3.
  50.     dipoleMomentLarge = 3.
  51.  
  52.     settings = {"autocorrelation":  True, \
  53.                 "densityAnalysis": False, \
  54.                 "vtf": False, \
  55.                 "dipoleAnalysis": False, \
  56.                 "radialDistribution": False,\
  57.                 "structureFactor": False }
  58.  
  59.  
  60.     # System parameters
  61.     nStepsEquil = 1000          #maximum amount of configs during equilibration
  62.     nConfigs = 10                #Amount of Configurations to simulate in equilibrium (for averaging measurements)
  63.     iStepsPerConfigEquil = 200  #integration steps between measurements during equilibration (energy)
  64.     iStepsPerConfigAna = 200    #integration steps between measurements during Analysis (exluding autocorrelation)
  65.     timeStep = 0.01
  66.     temperature = 1.
  67.     gamma = 1.
  68.  
  69.     #analysis parameters for multiple systems
  70.     nCycles = 20
  71.     dDensity = 0.04 #delta density | for densityAnalysis
  72.     dMoment = 4. #delta dipole moment | for dipoleAnalysis
  73.     iterateSmall = False
  74.     iterateLarge = True
  75.     if not settings["dipoleAnalysis"] and not settings["densityAnalysis"]: # only multiple runs if needed
  76.         nCycles = 1
  77.  
  78.  
  79.     #derived parameters / free parameters
  80.     boxLength = 60
  81.     boxVolume = np.power(boxLength,3)
  82.     volumeSmall = np.pi/6.*np.power(sigmaSmall,3)  #volume of small sphere
  83.     volumeLarge = np.pi/6.*np.power(sigmaLarge,3)  #volume of large sphere
  84.  
  85.     nSmall = int(densitySmall*boxVolume/volumeSmall)  #amount of small particles
  86.     nLarge = int(densityLarge*boxVolume/volumeLarge)  #amount of large particles
  87.  
  88.     wcaCut = 2.**(1. / 6.) # ~1.13
  89.  
  90.     particleTypes = []
  91.     small = ParticleType(type=0,n=nSmall,density=densitySmall,dipoleMoment=dipoleMomentSmall,sigma=sigmaSmall,epsilon=1,cut=wcaCut*sigmaSmall)
  92.     particleTypes.append(small)
  93.     large = ParticleType(type=1,n=nLarge,density=densityLarge,dipoleMoment=dipoleMomentLarge,sigma=sigmaLarge,epsilon=1,cut=wcaCut*sigmaLarge)
  94.     particleTypes.append(large)
  95.     nTypes = len(particleTypes)
  96.     nParticles = sum(particleType.n for particleType in particleTypes)
  97.  
  98.     #analysis Parameters
  99.     dc = DistanceCriterion(cut_off=particleTypes[1].cut) # cutoff big particle
  100.     cs = ClusterStructure(pair_criterion=dc)
  101.  
  102.     ################### Parameters end ###################
  103.  
  104.     if settings["densityAnalysis"]:
  105.         fpDensity = open('density{}.dat'.format(datetime.datetime.now()), mode="w+t")
  106.         fpDensity.write("densitySmall, densityLarge, Amount of clusters, Avg. cluster size\n")
  107.     if settings["dipoleAnalysis"]:
  108.         fpDipole = open('dipole{}.dat'.format(datetime.datetime.now()), mode="w+t")
  109.         fpDipole.write("dipoleMomentSmall, dipoleMomentLarge, Amount of clusters, Avg. cluster size\n")
  110.  
  111.     #init system #1
  112.     system = espressomd.System(box_l=boxLength*np.ones(3))
  113.     system.periodicity = [1, 1, 1] #periodic boundary condition in all directions
  114.     system.time_step = timeStep
  115.  
  116.     # Lennard-Jones interaction (WCA Potential)
  117.     for i in range(nTypes):
  118.         for j in range(i,nTypes):
  119.             lj_sig = combinationRuleSigma(particleTypes[i].sigma, particleTypes[j].sigma)
  120.             lj_cut = combinationRuleSigma(particleTypes[i].cut, particleTypes[j].cut)
  121.             lj_eps = combinationRuleEpsilon(particleTypes[i].epsilon, particleTypes[j].epsilon)
  122.             system.non_bonded_inter[particleTypes[i].type, particleTypes[j].type].lennard_jones.set_params(
  123.                                     epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift="auto")
  124.  
  125.     for iCycle in range(nCycles):
  126.         #init system #2
  127.         system.seed = iCycle
  128.         system.actors.clear()
  129.         #system.seed = 42
  130.  
  131.         #init particles
  132.         if system.part.n_part_types>0:
  133.             system.part.clear()
  134.         print("Particletypes:")
  135.         for p in particleTypes : print(p)
  136.         print("\n")
  137.         addParticles(particleTypes,system)
  138.         #print("nPart in system: ",system.number_of_particles(type=particleTypes[0].type))
  139.  
  140.         # Lennard Jones / WCA Equilibration - Cold, highly damped system
  141.         minDistS = 0
  142.         minDistL = 0
  143.         minDistSL = 0
  144.         cap = 1.0
  145.         system.force_cap=cap
  146.  
  147.         system.thermostat.turn_off()
  148.         system.thermostat.set_langevin(kT=temperature/10.0, gamma=gamma*5.0, seed=iCycle)
  149.         while minDistS< particleTypes[0].sigma or minDistL<particleTypes[1].sigma or minDistSL<(particleTypes[0].sigma+particleTypes[1].sigma)/2.0:
  150.             #Warmup Helper: Cap max. force, increase slowly for overlapping particles
  151.             minDistS = system.analysis.min_dist([particleTypes[0].type],[particleTypes[0].type])
  152.             minDistL = system.analysis.min_dist([particleTypes[1].type],[particleTypes[1].type])
  153.             minDistSL = system.analysis.min_dist([particleTypes[0].type],[particleTypes[1].type])
  154.             print(" - {:.2f}, {:.2f}, {:.2f}".format(minDistS,minDistL,minDistSL),end='\r')
  155.             cap += min(minDistS,minDistL,minDistSL)
  156.             system.force_cap=cap
  157.             system.integrator.run(10)
  158.  
  159.         system.force_cap=0 #turn off force cap
  160.         #print("adding magnetostatics to system...")
  161.         p3m = magnetostatics.DipolarP3M(prefactor=1, accuracy=1E-4)
  162.         system.actors.add(p3m)
  163.  
  164.  
  165.         ################### Energy Equilibration start ###################
  166.         print("Energy Equilibration ...")
  167.         system.thermostat.set_langevin(kT=temperature, gamma=gamma)
  168.         if settings["vtf"]:
  169.             fpVtf = open('trajectoryEquilibration_{}.vtf'.format(datetime.datetime.now()), mode='w+t')
  170.             _writevsf(system,fpVtf)
  171.         energy = system.analysis.energy()
  172.         energyOld = abs(energy['total'])
  173.         for i in range(1,nStepsEquil*iStepsPerConfigEquil):
  174.             if settings["vtf"]:
  175.                 _writevcf(system,fpVtf)
  176.             system.integrator.run(1)
  177.  
  178.             if i%iStepsPerConfigEquil==0:
  179.                 temp_measured = energy["kinetic"] / (1.5 * nParticles)
  180.                 energy = system.analysis.energy()
  181.                 energyNew = abs(energy['total'])
  182.                 relativeChange = abs(energyNew-energyOld)/energyOld
  183.                 print(" - t={0:.1f}, E={1:.2f},T={2:.4f},dE={3:.3f}".format(system.time, energyNew, temp_measured,relativeChange),end='\r')
  184.                 if relativeChange < 0.03 and energy['total']<0: #if energy change is small enough it's time for analysis
  185.                     print("\n - Old Energy {} New Energy {}".format(energyOld,energyNew))
  186.                     print(" - Temp Equilibration: break condition reached!")
  187.                     break
  188.                 energyOld=energyNew
  189.  
  190.         if settings["vtf"]:
  191.             fpVtf.close()
  192.         ################### Energy Equilibration end ###################
  193.  
  194.         if settings["densityAnalysis"]:
  195.             cs.run_for_all_pairs()
  196.             fpDensity.write("{0:f}, {1:f}, {2:d}, {3:f}\n".format(particleTypes[0].density, particleTypes[1].density, len(cs.clusters),averageClusterSize(cs.clusters)))
  197.             print(len(cs.clusters))
  198.         if settings["dipoleAnalysis"]:
  199.             cs.run_for_all_pairs()
  200.             fpDipole.write("{0:f}, {1:f}, {2:d}, {3:f}\n".format(particleTypes[0].dipoleMoment,particleTypes[1].dipoleMoment,len(cs.clusters),averageClusterSize(cs.clusters)))
  201.  
  202.  
  203.         # Integration
  204.         if settings["vtf"]:
  205.             fpVtf = open('trajectory_{}.vtf'.format(datetime.datetime.now()), mode='w+t')
  206.             _writevsf(system,fpVtf)
  207.         print("\nrunning simulation...")
  208.         system.time = 0.0
  209.         for i in range(nConfigs):
  210.             print(" - {} out of {} configurations done".format(i+1, nConfigs), end='\r')
  211.             system.integrator.run(iStepsPerConfigAna)
  212.             if settings["radialDistribution"]: #only averaging here
  213.                 system.analysis.append()
  214.             if settings["vtf"]:
  215.                 _writevcf(system,fpVtf)
  216.  
  217.         if settings["vtf"]:
  218.             fpVtf.close()
  219.  
  220.         system.time_step = timeStep
  221.         ################### Velocity autocorrelation start ###################
  222.         if settings["autocorrelation"]:
  223.             velObs = ParticleVelocities(ids=(0,))
  224.             cVel = Correlator(obs1=velObs, delta_N=1, corr_operation="scalar_product", tau_max=20.,tau_lin=16, compress1="discard1")
  225.             system.auto_update_accumulators.add(cVel)
  226.             system.integrator.run(1000)
  227.             cVel.finalize()
  228.             print(cVel.result())
  229.             #print(cVel.result().reshape([-1, 1]))
  230.             plt.ylim([0,2])
  231.             plt.plot(cVel.result()[:,0],cVel.result()[:,2],'bo')
  232.             plt.xlabel('$t$', fontsize=20)
  233.             plt.ylabel('$velocity autocorrelation$', fontsize=20)
  234.             plt.show()
  235.  
  236.         ################### Velocity autocorrelation end ###################
  237.  
  238.         ################### Radial distribution start ###################
  239.         if settings["radialDistribution"]:
  240.             bins = 100
  241.             rMin  = 0.0
  242.             rMax  = system.box_l[0]/4.0
  243.             r,rdfss = system.analysis.rdf(rdf_type='<rdf>',
  244.                                         type_list_a=[particleTypes[0].type],
  245.                                         type_list_b=[particleTypes[0].type],
  246.                                         r_min=rMin,
  247.                                         r_max=rMax,
  248.                                         r_bins=bins)
  249.             r,rdfsl = system.analysis.rdf(rdf_type='<rdf>',
  250.                                         type_list_a=[particleTypes[0].type],
  251.                                         type_list_b=[particleTypes[1].type],
  252.                                         r_min=rMin, r_max=rMax, r_bins=bins)
  253.             r,rdfll = system.analysis.rdf(rdf_type='<rdf>',
  254.                                         type_list_a=[particleTypes[1].type],
  255.                                         type_list_b=[particleTypes[1].type],
  256.                                         r_min=rMin, r_max=rMax, r_bins=bins)
  257.  
  258.             plt.figure(figsize=(10,6), dpi=80)
  259.             plt.plot(r[:],rdfss[:], label='Small')
  260.             plt.plot(r[:],rdfsl[:], label='SmallLarge')
  261.             plt.plot(r[:],rdfll[:], label='Large')
  262.             plt.xlabel('$r$', fontsize=20)
  263.             plt.ylabel('$g(r)$', fontsize=20)
  264.             plt.legend(fontsize=20)
  265.             plt.show()
  266.         ################### Radial distribution end ###################
  267.  
  268.         ################### Structure factor start ###################
  269.         if settings["structureFactor"]:
  270.             qs,ss = system.analysis.structure_factor(sf_types=[particleTypes[0].type],sf_order=60)
  271.             param, param_cov = curve_fit(func, qs, ss)
  272.             yss = (param[0]*np.exp(-param[1] * qs) + param[2]*(np.sin(param[3]*qs)) + param[4])
  273.  
  274.             ql,sl = system.analysis.structure_factor(sf_types=[particleTypes[1].type],sf_order=60)
  275.             param, param_cov = curve_fit(func, ql, sl)
  276.             ysl = (param[0]*np.exp(-param[1] * ql) + param[2]*(np.sin(param[3]*ql)) + param[4])
  277.  
  278.             qsl,ssl = system.analysis.structure_factor(sf_types=[particleTypes[0].type,particleTypes[1].type],sf_order=60)
  279.             param, param_cov = curve_fit(func, qsl, ssl)
  280.             yssl = (param[0]*np.exp(-param[1] * qsl) + param[2]*(np.sin(param[3]*qsl)) + param[4])
  281.  
  282.             fig, axs = plt.subplots(2, 2, sharex=True, sharey=True)
  283.             axs[0, 0].plot(qs[:],ss[:],'bx', alpha=0.5)
  284.             #plt.plot(qs[:],ss[:],'bx', label='Small', alpha=0.5)
  285.             axs[0, 0].plot(qs[:], yss[:], 'b-')
  286.             axs[0, 0].set_title('Small')
  287.  
  288.             axs[0, 1].plot(ql[:],sl[:],'gx', alpha=0.5)
  289.             axs[0, 1].plot(ql[:], ysl[:], 'g-')
  290.             axs[0, 1].set_title('Large')
  291.  
  292.  
  293.             axs[1, 0].plot(qsl[:],ssl[:],'rx', alpha=0.5)
  294.             axs[1, 0].plot(qsl[:], yssl[:], 'r-')
  295.             axs[1, 0].set_title('Small & Large')
  296.  
  297.             axs[1, 1].plot(qs[:], yss[:], 'b-', label='S fit')
  298.             axs[1, 1].plot(ql[:], ysl[:], 'g-', label='L fit')
  299.             axs[1, 1].plot(qsl[:], yssl[:], 'r-', label='S&L fit')
  300.             axs[1, 1].legend(loc="upper right")
  301.             axs[1, 1].set_title('All fits')
  302.  
  303.             for ax in axs.flat:
  304.                 ax.set(xlabel='$q$', ylabel='$S(q)$')
  305.  
  306.             for ax in axs.flat:
  307.                 ax.label_outer()
  308.  
  309.             #axs[1,1].
  310.             #axs[1,1].
  311.             plt.show()
  312.         ################### Structure factor end ###################
  313.  
  314.  
  315.         if settings["dipoleAnalysis"]:
  316.             if(iterateSmall):
  317.                 particleTypes[0].dipoleMoment += dMoment/nCycles
  318.             if(iterateLarge):
  319.                 particleTypes[1].dipoleMoment += dMoment/nCycles
  320.  
  321.         if settings["densityAnalysis"]:
  322.             if(iterateSmall):
  323.                 particleTypes[0].density += dDensity/nCycles
  324.                 particleTypes[0].n = int(particleTypes[0].density*boxVolume/volumeSmall)
  325.             if(iterateLarge):
  326.                 particleTypes[1].density += dDensity/nCycles
  327.                 particleTypes[1].n = int(particleTypes[1].density*boxVolume/volumeLarge)
  328.  
  329.  
  330.     if settings["densityAnalysis"]:
  331.         fpDensity.write("\n\n\nSystem Parameters:\nTemperature: {0:f}\n".format(temperature))
  332.         for p in particleTypes:
  333.              fpDensity.write(str(p)+'\n')
  334.         fpDensity.close()
  335.     elif settings["dipoleAnalysis"]:
  336.         fpDipole.write("\n\n\nSystem Parameters:\n")
  337.         for p in particleTypes:
  338.              fpDipole.write(str(p)+'\n')
  339.         fpDipole.close()
  340.  
  341.  
  342.  
  343.     # Analysis
  344.     print("\n\nANALYSIS:")
  345.     cs.run_for_all_pairs()
  346.     print("  #Clusters: {:d}".format(len(cs.clusters)))
  347.     max_size = 0;
  348.     nSmallInCluster = 0
  349.     nLargeInCluster = 0
  350.     for c in cs.clusters:
  351.         print("  cluster id: {}".format(c[0]))
  352.         nSmall = 0
  353.         nLarge = 0
  354.         for p in c[1].particles():
  355.             if p.type == particleTypes[0].type:
  356.                 nSmall+=1
  357.             elif p.type == particleTypes[1].type:
  358.                 nLarge+=1
  359.         nSmallInCluster+=nSmall
  360.         nLargeInCluster+=nLarge
  361.         print("    nSmall in cluster: {} \n    nLarge in cluster: {}".format(nSmall,nLarge))
  362.  
  363.         if max_size<c[1].size():
  364.             max_size = c[1].size()
  365.         fd = c[1].fractal_dimension(dr=0.1)[0]
  366.         if fd == fd: # if fd is not nan
  367.             print("    Fractal dimension:",fd)
  368.     print("  biggest cluster: {:d}".format(max_size))
  369.     print("  {}/{} small particles in cluster".format(nSmallInCluster,particleTypes[0].n))
  370.     print("  {}/{} large particles in cluster".format(nLargeInCluster,particleTypes[1].n))
  371.  
  372.  
  373.  
  374. if __name__ == "__main__":
  375.     main()
  376.  
  377.  
  378. #backup code
  379. #system.time_step = timeStep
  380. ################### Velocity autocorrelation start ###################
  381. #velObs = espressomd.observables.ParticleVelocities(ids=(0,))
  382. #accumulator = espressomd.accumulators.TimeSeries(obs=velObs, delta_N=1)
  383. #system.auto_update_accumulators.add(accumulator)
  384. #system.integrator.run(300000)
  385. #print(accumulator.time_series())
  386. #acf = tidynamics.acf(accumulator.time_series())
  387. #print(acf)
  388. #fig = plt.figure()
  389. #ax = fig.gca()
  390. #ax.set_xticks(np.arange(0, 1000, 100))
  391. #ax.set_yticks(np.arange(min(acf), max(acf), 0.5))
  392. #plt.plot(acf)
  393. #plt.grid()
  394. #plt.show()
  395. #plt.plot(acf) # Plot list. x-values assumed to be [0, 1, 2, 3]
  396. #plt.show()  # Optional when using IPython interpreter
  397. #vector_autocorrelate(accumulator.time_series())
  398. ################### Velocity autocorrelation end ###################
  399.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!