In [1]:

    

%load_ext autoreload
%autoreload 1


    

In [2]:

    

%aimport pyparm
%aimport pyparm.packmin
%aimport parview
import pyparm.d3 as sim
from pyparm.util import norm
from math import pi
from itertools import count
import jsonpickle


    

In [3]:

    

N = 32
alpha = 0.5
L0 = 40.

#diameter is 1
V_cap = pi / 6.
V_cyl = alpha * pi / 4.
V_SC = N * (V_cap + V_cyl)
V = L0**3
print(V_SC, V_SC / V)


    

    




29.321531433504735 0.00045814892864851146

In [4]:

    

scavec = sim.SCatomvec([1.0] * N)
springs = sim.SCSpringList(scavec, 1.0, 1.0, alpha)
box = sim.OriginBox(L0)
for a in scavec:
    a.x = box.randLoc()
    a.v = sim.vec()
    a.f = sim.vec()

constraints = []
for a1, a2 in scavec.all_pairs():
    dc = sim.distConstraint(a1, a2, alpha)
    dc.apply_positions(box)
    constraints.append(dc)
print(springs.energy(box))


    

    




0.0

In [5]:

    

collec = sim.collectionNLCG(box, scavec, 0.01, 1e-4, [springs], [], constraints)


    

In [6]:

    

m = 1000
locs = []
Ls = []

#while True:
for i in range(m):
    for _ in range(200):
        collec.timestep()
    fdotf = collec.fdotf() / N
    if i % 10 == 0:
        H = collec.Hamiltonian()
        V = box.V()
        phi = V_SC / V
        print('{:4d}: {:10.4g}, phi {:12.9f}. H={:15.12f}  (L: {:9.3g} / {:9.3g})'.format(
                i, fdotf, phi, H, box.L(), ((1 + alpha) * 2)))
    locs.append([a.x for a in scavec])
    Ls.append(box.L())
    if fdotf < 1e-12:
        print('Done.')
        break
    if box.L() <= (1 + alpha) * 2:
        print('Box shrunk too much.')
        break


    

    




   0:   0.008129, phi  0.000574915. H= 5.100148778198  (L:      37.1 /         3)
  10:  0.0008792, phi  0.001748103. H= 1.677334681017  (L:      25.6 /         3)
  20:  0.0003143, phi  0.002919834. H= 1.004223808740  (L:      21.6 /         3)
  30:  0.0001601, phi  0.004091214. H= 0.716697115584  (L:      19.3 /         3)
  40:  9.659e-05, phi  0.005261629. H= 0.557273565356  (L:      17.7 /         3)
  50:   4.61e-05, phi  0.007603219. H= 0.385648489597  (L:      15.7 /         3)
  60:   2.15e-05, phi  0.011103051. H= 0.264087504756  (L:      13.8 /         3)
  70:  1.242e-05, phi  0.014596504. H= 0.200882010192  (L:      12.6 /         3)
  80:  7.393e-06, phi  0.018922700. H= 0.154955162393  (L:      11.6 /         3)
  90:  3.603e-06, phi  0.027068778. H= 0.108322863296  (L:      10.3 /         3)
 100:  2.102e-06, phi  0.035182225. H= 0.083342624525  (L:      9.41 /         3)
 110:  1.197e-06, phi  0.046577630. H= 0.062952370310  (L:      8.57 /         3)
 120:  6.263e-07, phi  0.063748146. H= 0.045996306199  (L:      7.72 /         3)
 130:  3.787e-07, phi  0.080726678. H= 0.036322466065  (L:      7.13 /         3)
 140:  1.926e-07, phi  0.111319686. H= 0.026340274304  (L:      6.41 /         3)
 150:  1.101e-07, phi  0.145180877. H= 0.020196807958  (L:      5.87 /         3)
 160:  6.495e-08, phi  0.186759843. H= 0.015700687234  (L:      5.39 /         3)
 170:  3.134e-08, phi  0.258334141. H= 0.011350830811  (L:      4.84 /         3)
 180:  1.781e-08, phi  0.332101578. H= 0.008829734359  (L:      4.45 /         3)
 190:  1.106e-08, phi  0.436963827. H= 0.006711866764  (L:      4.06 /         3)
 200:  7.473e-09, phi  0.557509754. H= 0.005261548753  (L:      3.75 /         3)
 210:  7.352e-10, phi  0.668411508. H= 0.004387760934  (L:      3.53 /         3)
 220:  1.798e-10, phi  0.681041843. H= 0.004306487240  (L:       3.5 /         3)
 230:  6.365e-12, phi  0.686143988. H= 0.004274311728  (L:       3.5 /         3)
 240:  3.419e-11, phi  0.686422192. H= 0.004272582175  (L:       3.5 /         3)
 250:   1.48e-12, phi  0.689195384. H= 0.004255327797  (L:      3.49 /         3)
Done.

In [7]:

    

locs = np.asarray(locs)
Ls = np.asarray(Ls)
print(locs.shape, Ls.shape)


    

    




(259, 64, 3) (259,)

In [39]:

    

dx = 0.5
x2 = np.sqrt((4 - dx**2)/2.)

np.sqrt(x2**2 + x2**2 + dx**2)


    

    
Out[39]:





2.0

In [83]:

    

locs1, locs2 = locs[:, ::2, :], locs[:, 1::2, :]
rs = (locs1 + locs2) / 2.
axes = (locs2 - locs1)
names = [[str(i+1)] for i in range(N)]
diameters = [1.0]*N

frames = []
for n, L, locset, rset, axset in zip(count(), Ls, locs, rs, axes):
    scs = [parview.Spherocylinder(r/L % 1.0 - 0.5, d/L, ax/L * (alpha + 1.) / alpha, nm) 
           for r, d, ax, nm in zip(rset, diameters, axset, names)]
    
    text='%d, phi=%5.3f' % (n, V_SC / (L**3))
    frames.append(parview.Frame(spherocylinders=scs, spheres=sphs, text=text))


    

In [85]:

    

jsontext = jsonpickle.encode(frames, unpicklable=False)
with open('spherocylinders.json', 'w') as f:
    f.write(jsontext)


    

In [86]:

    

collec.potentialenergy()


    

    
Out[86]:





8.731006552961472e-07

In [ ]:

    

for