Commit 9fca17f4 authored by Robin Worreby's avatar Robin Worreby
Browse files

Merge remote-tracking branch 'upstream/master'

parents 82d599c9 8dab325d
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
import matplotlib.pyplot as plt
import numpy as np
from scipy.constants import physical_constants
import time
def unique(l):
newk = []
for i in l:
if i not in newk:
newk.append(i)
return newk
### periodic boundary conditions
def pbc(v,nx,ny):
b=np.array(v)
dim=len(b.shape)
if dim==2:
for a in v:
if a[0] <0:
a[0]=a[0]+nx
elif a[0] > nx -1 :
a[0]=a[0]-nx
if a[1] <0:
a[1]=a[1]+ny
elif a[1] > ny -1:
a[1]=a[1]-ny
else:
if v[0] <0:
v[0]=v[0]+nx
elif v[0] > nx -1 :
v[0]=v[0]-nx
if v[1] <0:
v[1]=v[1]+ny
elif v[1] > ny -1:
v[1]=v[1]-ny
return v
### map lattice system coordinate to cartesian coordinates
def coordinates(m,dx,dy):
x=[]
y=[]
c=[]
for i in m:
x.append(i[0]*dx+i[2]*dx/2)
y.append(i[1]*dy-i[2]*dy/2)
if i[3]==0:
c.append('blue')
else:
c.append('red')
return x,y,c
### list of possible events with reference to the lattice coordinate convention
### the forth index means unbound (0) or bound (1) molecule
def events(m,selection,nx,ny):
allevents=[]
set=[]
#### list of first neighbors (relative position)
# e.g.
# set[0] == (6)neigh. at site 0 (ids==0)
# set[1] == (6)neigh. at site 1 (ids==1)
set.append([[0,0,1,0],[1,0,0,0],[0,1,1,0],[-1,1,1,0],[-1,0,0,0],[-1,0,1,0]])
set.append([[1,-1,-1,0],[1,0,0,0],[1,0,-1,0],[0,0,-1,0],[-1,0,0,0],[0,-1,-1,0]])
for i in selection:
#### le is the list of possible events
le=[]
#### consider only molecules that are not binded
if m[i][3] == 0:
#### the funcion neighbors returns "n" the number of occupied
#### 1st nb sites and "l" implicit and explicit list of
#### the 1st nb sites occupied e.g. (2,[[0,[3,3,0,0]],[5,[2,3,1,0]]])
n,l=neighbors(m,i,nx,ny)
ids=m[i][2]
#### if there are no first neigh then 6 diffusions possible
#### add all possible displacements (tobeadded) for molecule i (involved)
thepoint=np.array(m[i])
if n==0:
tobeadded=pbc((thepoint + np.array(set[int(ids)])).tolist(),nx,ny)
involved=[i]
for d in tobeadded:
le.append([[d],'d',involved])
#le=le+pbc(tobeadded,nx,ny)
#le=[c+['d']+involved for c in le]
else:
# n = ... 1 |2 |3 |4 |5 ,where |=='or'
away=[0,1,2,3,4,5]
dr=[]
dstring='d'+str(n)
if n==2:
ndist=np.abs(l[0][0]-l[1][0])
if ndist==1 or ndist==5:
drstring='drl'+str(n) # 2 neighbooring nb
else:
drstring='dr'+str(n) # 2 non-neighbooring nb
else:
drstring='dr'+str(n)
# n = ... 1 |3 |4 |5 ,where |=='or'
for nn in l:
ni=nn[0]
ns=nn[1]
away.remove(ni)
#### binding
involved=[i,m.index(ns)]
new1=list(m[i])
new1[3]=1
new2=list(ns)
new2[3]=1
le.append([[new1,new2],'b',involved])
#### diffusion around the neighbor
for aw in [1,5]:
nsite=np.mod(ni+aw,6)
if nsite not in dr:
dr.append(nsite)
for nn in l:
ni=nn[0]
#dr = list(filter((ni).__ne__, dr))
if ni in dr:
dr.remove(ni)
#### for all possible mooves (dr) ![remaining around at least one neighbor] append event.
#### The event is coded as 'dri' where 'i' is the position relative to m[i]
for nsite in dr:
if nsite in away:
away.remove(nsite)
aux=(thepoint + np.array(set[int(ids)][nsite])).tolist()
aux=pbc(aux,nx,ny)
le.append([[aux],drstring,[i]])
#### away contains the possible moves to sites non-neighbooring to any nb
#### diffusion away
for nsite in away:
dstring='d'+str(n)
aux=(thepoint + np.array(set[int(ids)][nsite])).tolist()
aux=pbc(aux,nx,ny)
le.append([[aux],dstring,[i]])
allevents=allevents+le
return allevents
### compute the total rate of events
def total_rate(events,rates):
R=np.sum([rates[i[1]] for i in events])
return R
### select randomly an event from the list of possible events
def find_event(R,rates,events):
sum__l=0
rho1=np.random.random()
target=rho1*R
sum_l=rates[events[0][1]]
if target<sum_l:
the_event=events[0]
else:
for i in events[1:]:
#sum_U=sum_l+rates[i[1]]
sum_l+=rates[i[1]]
if target<sum_l:
the_event=i
break
#sum_l=sum_U
return the_event
### neighbors of a molecule according to the lattice coordinate convention
def neighbors(a,id,nx,ny):
nneig=0
allneig=[]
f=[]
f.append([[0,0,1,0],[1,0,0,0],[0,1,1,0],[-1,1,1,0],[-1,0,0,0],[-1,0,1,0]])
f.append([[1,-1,-1,0],[1,0,0,0],[1,0,-1,0],[0,0,-1,0],[-1,0,0,0],[0,-1,-1,0]])
thepoint=np.array(a[id])
ids=a[id][2]
for n in range(6):
theneig=(thepoint + np.array(f[int(ids)][n])).tolist()
theneig=pbc(theneig,nx,ny)
alt1=list(theneig)
alt2=list(theneig)
alt1[3]=1
alt2[3]=0
if alt1 in a:
nneig=nneig+1
allneig.append([n,alt1])
if alt2 in a:
nneig=nneig+1
allneig.append([n,alt2])
return nneig,allneig
### function to execute an event
def apply_event(molecules,selected_event,possible_events,nx,ny):
new_positions=selected_event[0]
tobeupdated=[]
theinvolved=0
id_involved=selected_event[2]
for involved in id_involved:
tobeupdated.append(involved)
nu,lu=neighbors(molecules,involved,nx,ny)
for ii in lu:
theindex=molecules.index(ii[1])
if theindex not in tobeupdated:
tobeupdated.append(theindex)
for involved in id_involved:
molecules[involved]=new_positions[theinvolved]
theinvolved=theinvolved+1
theinvolved=0
for involved in id_involved:
nu,lu=neighbors(molecules,involved,nx,ny)
for ii in lu:
theindex=molecules.index(ii[1])
if theindex not in tobeupdated:
tobeupdated.append(theindex)
theinvolved=theinvolved+1
i_shift = 0
oldevents=list(possible_events)
for i_oe, oe in enumerate(possible_events):
for tu in tobeupdated:
if tu in oe[2]:
del oldevents[i_oe+i_shift]
i_shift -= 1
break
possible_events=list(oldevents)
# oldevents=list(possible_events)
# for oe in possible_events:
# for tu in tobeupdated:
# if tu in oe[2]:
# oldevents.remove(oe)
# break
# possible_events=list(oldevents)
return possible_events,tobeupdated
\ No newline at end of file
import matplotlib.pyplot as plt
import numpy as np
from scipy.constants import physical_constants
#### this subroutine identifies all molecules connected to a given one
#### simple BUT NOT EFFICIENT
def allconnected(m,id,nx,ny):
theconnected=[m[id]]
n,l=neighbors(m,id,nx,ny)
if n>0:
remain=list(m)
remain.remove(m[id])
for i in l:
#theconnected.append(m[i])
remain.remove(m[i])
for i in l:
remain.append(m[i])
theconnected=theconnected+allconnected(remain,len(remain)-1,nx,ny)
for j in theconnected:
if j in remain:
remain.remove(j)
return theconnected
def ene(mol,de,nx,ny):
totene=0.0
for m in range(len(mol)):
n,l =neighbors(mol,m,nx,ny)
totene=totene + n*de
return totene/2
#### PERIODIC BOUNDARY CONTITIONS
def pbc(a,nx,ny):
if a[0] <0:
a[0]=a[0]+nx
elif a[0] > nx -1 :
a[0]=a[0]-nx
if a[1] <0:
a[1]=a[1]+ny
elif a[1] > ny -1:
a[1]=a[1]-ny
return a
#### CONVERT the lattice coordinates in cartesian coordinates
#### used for plotting
def coordinates(m,dx,dy):
x=[]
y=[]
for i in m:
x.append(i[0]*dx+i[2]*dx/2)
y.append(i[1]*dy-i[2]*dy/2)
return x,y
#### find the 1st neighboring molecules to a given one
def neighbors(a,id,nx,ny):
nneig=0
allneig=[]
thepoint=np.array(a[id])
if a[id][2] == 0:
for inc in [[0,0,1],[1,0,0],[0,1,1],[-1,1,1],[-1,0,0],[-1,0,1]]:
theneig=(thepoint + np.array(inc)).tolist()
theneig=pbc(theneig,nx,ny)
if theneig in a:
nneig=nneig+1
allneig.append(a.index(theneig))
else:
for inc in [[1,-1,-1],[1,0,0],[1,0,-1],[0,0,-1],[-1,0,0],[0,-1,-1]]:
theneig=(thepoint + np.array(inc)).tolist()
theneig=pbc(theneig,nx,ny)
if theneig in a:
nneig=nneig+1
allneig.append(a.index(theneig))
return nneig,allneig
\ No newline at end of file
This diff is collapsed.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment