-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathInitconfig.f90
156 lines (145 loc) · 4.97 KB
/
Initconfig.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
SUBROUTINE Initconfig(x1,y1,z1,boa,coa,aoc,boc,Rs1,Rodl,N,Nc)
!!This Subroutine initializes the position of beads with length of rod scaled to ellipsoidal (cavity) radius
!!Rod length (Cavity) radius is chosen to be the length scale in the problem
IMPLICIT NONE
INTEGER*8 :: i,N,Nc,j
REAL*8 :: px,py,pz,pnorm,Rs,Rs1,dr,xnorm,x1(N,Nc),y1(N,Nc),z1(N,Nc),Rodl,boa,coa,aoc,boc,&
R_cx(Nc),R_cy(Nc),R_cz(Nc),num,dr1
REAL*4 :: rand
dr = Rs1*0.05d0
num = 3.d0
R_cx(1) = -Rs1/num
R_cy(1) = 0.d0
R_cz(1) = Rs1/num
R_cx(2) = Rs1/num
R_cy(2) = 0.d0
R_cz(2) = -Rs1/num
R_cx(3) = Rs1/num
R_cy(3) = 0.d0
R_cz(3) = Rs1/num
R_cx(4) = -Rs1/num
R_cy(4) = 0.d0
R_cz(4) = -Rs1/num
Rs = Rs1!/2.d0
j = 1
xnorm = 1000.d0
x1(1:N,j) = 1.d0
y1(1:N,j) = 1000.d0
z1(1:N,j) = -1.d0
dr1 = 0.8d0
DO WHILE(xnorm .GT. (Rs - dr) .OR. x1(1,j) .GT. -dr1 .OR. z1(1,j) .LT. dr1)
x1(1,j) = Rs*(rand(0) - 0.5d0)
y1(1,j) = Rs*(rand(0) - 0.5d0)
z1(1,j) = Rs*(rand(0) - 0.5d0)
xnorm = DSQRT((x1(1,j)/aoc-R_cx(j))**2.d0+(y1(1,j)/boc-R_cy(j))**2.d0+(z1(1,j)-R_cz(j))**2.d0)
ENDDO
DO i = 1,N-1
xnorm = 1000.d0
DO WHILE(xnorm .GT. (Rs - dr) .OR. x1(i+1,j) .GT. -dr1 .OR. z1(i+1,j) .LT. dr1)
px = (rand(0) - 0.5d0)*2.d0
py = (rand(0) - 0.5d0)*2.d0
pz = (rand(0) - 0.5d0)*2.d0
pnorm = DSQRT(px**2.d0 + py**2.d0 + pz**2.d0)
px = Rodl*px/pnorm
py = Rodl*py/pnorm
pz = Rodl*pz/pnorm
x1(i+1,j) = x1(i,j) + px
y1(i+1,j) = y1(i,j) + py
z1(i+1,j) = z1(i,j) + pz
xnorm = DSQRT((x1(i+1,j)/aoc)**2.d0+(y1(i+1,j)/boc)**2.d0+(z1(i+1,j))**2.d0)
!PRINT*,i,j,x1(i+1,j),y1(i+1,j),z1(i+1,j),xnorm
ENDDO
!PRINT*,i+1,j,x1(i+1,j),y1(i+1,j),z1(i+1,j)
!PAUSE
ENDDO
j = 2
xnorm = 1000.d0
x1(1:N,j) = -1.d0
y1(1:N,j) = 1000.d0
z1(1:N,j) = 1.d0
DO WHILE(xnorm .GT. (Rs - dr) .OR. x1(1,j) .LT. dr1 .OR. z1(1,j) .GT. -dr1)
x1(1,j) = Rs*(rand(0) - 0.5d0)
y1(1,j) = Rs*(rand(0) - 0.5d0)
z1(1,j) = Rs*(rand(0) - 0.5d0)
xnorm = DSQRT((x1(1,j)/aoc-R_cx(j))**2.d0+(y1(1,j)/boc-R_cy(j))**2.d0+(z1(1,j)-R_cz(j))**2.d0)
ENDDO
DO i = 1,N-1
xnorm = 1000.d0
DO WHILE(xnorm .GT. (Rs - dr) .OR. x1(i+1,j) .LT. dr1 .OR. z1(i+1,j) .GT. -dr1)
px = (rand(0) - 0.5d0)*2.d0
py = (rand(0) - 0.5d0)*2.d0
pz = (rand(0) - 0.5d0)*2.d0
pnorm = DSQRT(px**2.d0 + py**2.d0 + pz**2.d0)
px = Rodl*px/pnorm
py = Rodl*py/pnorm
pz = Rodl*pz/pnorm
x1(i+1,j) = x1(i,j) + px
y1(i+1,j) = y1(i,j) + py
z1(i+1,j) = z1(i,j) + pz
xnorm = DSQRT((x1(i+1,j)/aoc)**2.d0+(y1(i+1,j)/boc)**2.d0+(z1(i+1,j))**2.d0)
!PRINT*,i,j,x1(i+1,j),y1(i+1,j),z1(i+1,j),xnorm
ENDDO
!PRINT*,i+1,j,x1(i+1,j),y1(i+1,j),z1(i+1,j)
!PAUSE
ENDDO
j = 3
xnorm = 1000.d0
x1(1:N,j) = -1.d0
y1(1:N,j) = 1000.d0
z1(1:N,j) = -1.d0
DO WHILE(xnorm .GT. (Rs - dr) .OR. x1(1,j) .LT. dr1 .OR. z1(1,j) .LT. dr1)
x1(1,j) = Rs*(rand(0) - 0.5d0)
y1(1,j) = Rs*(rand(0) - 0.5d0)
z1(1,j) = Rs*(rand(0) - 0.5d0)
xnorm = DSQRT((x1(1,j)/aoc-R_cx(j))**2.d0+(y1(1,j)/boc-R_cy(j))**2.d0+(z1(1,j)-R_cz(j))**2.d0)
ENDDO
DO i = 1,N-1
xnorm = 1000.d0
DO WHILE(xnorm .GT. (Rs - dr) .OR. x1(i+1,j) .LT. dr1 .OR. z1(i+1,j) .LT. dr1)
px = (rand(0) - 0.5d0)*2.d0
py = (rand(0) - 0.5d0)*2.d0
pz = (rand(0) - 0.5d0)*2.d0
pnorm = DSQRT(px**2.d0 + py**2.d0 + pz**2.d0)
px = Rodl*px/pnorm
py = Rodl*py/pnorm
pz = Rodl*pz/pnorm
x1(i+1,j) = x1(i,j) + px
y1(i+1,j) = y1(i,j) + py
z1(i+1,j) = z1(i,j) + pz
xnorm = DSQRT((x1(i+1,j)/aoc)**2.d0+(y1(i+1,j)/boc)**2.d0+(z1(i+1,j))**2.d0)
!PRINT*,i,j,x1(i+1,j),y1(i+1,j),z1(i+1,j),xnorm
ENDDO
!PRINT*,i+1,j,x1(i+1,j),y1(i+1,j),z1(i+1,j)
!PAUSE
ENDDO
j = 4
xnorm = 1000.d0
x1(1:N,j) = 1.d0
y1(1:N,j) = 1000.d0
z1(1:N,j) = 1.d0
DO WHILE(xnorm .GT. (Rs - dr) .OR. x1(1,j) .GT. -dr1 .OR. z1(1,j) .GT. -dr1)
x1(1,j) = Rs*(rand(0) - 0.5d0)
y1(1,j) = Rs*(rand(0) - 0.5d0)
z1(1,j) = Rs*(rand(0) - 0.5d0)
xnorm = DSQRT((x1(1,j)/aoc-R_cx(j))**2.d0+(y1(1,j)/boc-R_cy(j))**2.d0+(z1(1,j)-R_cz(j))**2.d0)
ENDDO
DO i = 1,N-1
xnorm = 1000.d0
DO WHILE(xnorm .GT. (Rs - dr) .OR. x1(i+1,j) .GT. -dr1 .OR. z1(i+1,j) .GT. -dr1)
px = (rand(0) - 0.5d0)*2.d0
py = (rand(0) - 0.5d0)*2.d0
pz = (rand(0) - 0.5d0)*2.d0
pnorm = DSQRT(px**2.d0 + py**2.d0 + pz**2.d0)
px = Rodl*px/pnorm
py = Rodl*py/pnorm
pz = Rodl*pz/pnorm
x1(i+1,j) = x1(i,j) + px
y1(i+1,j) = y1(i,j) + py
z1(i+1,j) = z1(i,j) + pz
xnorm = DSQRT((x1(i+1,j)/aoc)**2.d0+(y1(i+1,j)/boc)**2.d0+(z1(i+1,j))**2.d0)
!PRINT*,i,j,x1(i+1,j),y1(i+1,j),z1(i+1,j),xnorm
ENDDO
!PRINT*,i+1,j,x1(i+1,j),y1(i+1,j),z1(i+1,j)
!PAUSE
ENDDO
END SUBROUTINE Initconfig