CasualGamer

COPHY_Main

Jan 10th, 2021
953
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  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.  
RAW Paste Data