ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/calc_shapes.F90
(Generate patch)

Comparing trunk/OOPSE/libmdtools/calc_shapes.F90 (file contents):
Revision 1320 by gezelter, Tue Jun 29 21:15:22 2004 UTC vs.
Revision 1321 by gezelter, Fri Jul 2 21:41:10 2004 UTC

# Line 24 | Line 24 | module shapes
24  
25    public :: do_shape_pair
26  
27 <  type :: ShapeList
28 <     integer :: nContactFuncs = 0
29 <     integer :: nRangeFuncs = 0
30 <     integer :: nStrengthFuncs = 0
31 <     integer :: bigL = 0
32 <     integer :: bigM = 0
33 <     integer, allocatable, dimension(:) :: ContactFuncLValue
34 <     integer, allocatable, dimension(:) :: ContactFuncMValue
35 <     integer, allocatable, dimension(:) :: ContactFunctionType
36 <     real(kind=dp), allocatable, dimension(:) :: ContactFuncCoefficient
37 <     integer, allocatable, dimension(:) :: RangeFuncLValue
38 <     integer, allocatable, dimension(:) :: RangeFuncMValue
39 <     integer, allocatable, dimension(:) :: RangeFunctionType
40 <     real(kind=dp), allocatable, dimension(:) :: RangeFuncCoefficient
41 <     integer, allocatable, dimension(:) :: StrengthFuncLValue
42 <     integer, allocatable, dimension(:) :: StrengthFuncMValue
43 <     integer, allocatable, dimension(:) :: StrengthFunctionType
44 <     real(kind=dp), allocatable, dimension(:) :: StrengthFuncCoefficient
45 <     logical :: isLJ = .false.
46 <     real ( kind = dp )  :: epsilon = 0.0_dp
47 <     real ( kind = dp )  :: sigma = 0.0_dp        
27 >
28 >  type, private :: Shape
29 >     integer :: atid
30 >     integer :: nContactFuncs
31 >     integer :: nRangeFuncs
32 >     integer :: nStrengthFuncs
33 >     integer :: bigL
34 >     integer :: bigM
35 >     integer, pointer, dimension(:) :: ContactFuncLValue             => null()
36 >     integer, pointer, dimension(:) :: ContactFuncMValue             => null()
37 >     integer, pointer, dimension(:) :: ContactFunctionType           => null()
38 >     real(kind=dp), pointer, dimension(:) :: ContactFuncCoefficient  => null()
39 >     integer, pointer, dimension(:) :: RangeFuncLValue               => null()
40 >     integer, pointer, dimension(:) :: RangeFuncMValue               => null()
41 >     integer, pointer, dimension(:) :: RangeFunctionType             => null()
42 >     real(kind=dp), pointer, dimension(:) :: RangeFuncCoefficient    => null()
43 >     integer, pointer, dimension(:) :: StrengthFuncLValue            => null()
44 >     integer, pointer, dimension(:) :: StrengthFuncMValue            => null()
45 >     integer, pointer, dimension(:) :: StrengthFunctionType          => null()
46 >     real(kind=dp), pointer, dimension(:) :: StrengthFuncCoefficient => null()
47 >     logical :: isLJ
48 >     real ( kind = dp )  :: epsilon
49 >     real ( kind = dp )  :: sigma
50 >  end type Shape
51 >  
52 >  type, private :: ShapeList
53 >     integer :: n_shapes = 0
54 >     integer :: currentShape = 0
55 >     type (Shape), pointer :: Shapes(:)      => null()
56 >     integer, pointer      :: atidToShape(:) => null()
57    end type ShapeList
58 +  
59 +  type(ShapeList), save :: ShapeMap
60  
50 –  type(ShapeList), dimension(:),allocatable :: ShapeMap
51 –
61    integer :: lmax
62    real (kind=dp), allocatable, dimension(:,:) :: plm_i, dlm_i, plm_j, dlm_j
63    real (kind=dp), allocatable, dimension(:) :: tm_i, dtm_i, um_i, dum_i
# Line 56 | Line 65 | contains  
65  
66   contains  
67  
68 <  subroutine createShapeMap(status)
69 <    integer :: nAtypes
68 >  subroutine newShapeType(nContactFuncs, ContactFuncLValue, &
69 >       ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
70 >       nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
71 >       RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
72 >       StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
73 >       myAtid, status)
74 >
75 >    integer :: nContactFuncs
76 >    integer :: nRangeFuncs
77 >    integer :: nStrengthFuncs
78 >    integer :: shape_ident
79      integer :: status
80 <    integer :: i
81 <    real (kind=DP) :: thisDP
82 <    logical :: thisProperty
80 >    integer :: myAtid
81 >    integer :: bigL
82 >    integer :: bigM
83 >    integer :: j, me, nShapeTypes, nLJTypes, ntypes, current, alloc_stat
84 >    integer, pointer :: MatchList(:) => null()
85  
86 +    integer, dimension(nContactFuncs) :: ContactFuncLValue          
87 +    integer, dimension(nContactFuncs) :: ContactFuncMValue          
88 +    integer, dimension(nContactFuncs) :: ContactFunctionType        
89 +    real(kind=dp), dimension(nContactFuncs) :: ContactFuncCoefficient
90 +    integer, dimension(nRangeFuncs) :: RangeFuncLValue            
91 +    integer, dimension(nRangeFuncs) :: RangeFuncMValue            
92 +    integer, dimension(nRangeFuncs) :: RangeFunctionType          
93 +    real(kind=dp), dimension(nRangeFuncs) :: RangeFuncCoefficient  
94 +    integer, dimension(nStrengthFuncs) :: StrengthFuncLValue          
95 +    integer, dimension(nStrengthFuncs) :: StrengthFuncMValue          
96 +    integer, dimension(nStrengthFuncs) :: StrengthFunctionType        
97 +    real(kind=dp), dimension(nStrengthFuncs) :: StrengthFuncCoefficient
98 +
99      status = 0
100 +    ! check to see if this is the first time into this routine...
101 +    if (.not.associated(ShapeMap%Shapes)) then
102  
103 <    nAtypes = getSize(atypes)
103 >       call getMatchingElementList(atypes, "is_Shape", .true., &
104 >            nShapeTypes, MatchList)
105 >      
106 >       call getMatchingElementList(atypes, "is_LJ", .true., &
107 >            nLJTypes, MatchList)
108 >      
109 >       ShapeMap%n_shapes = nShapeTypes + nLJTypes
110 >      
111 >       allocate(ShapeMap%Shapes(nShapeTypes + nLJTypes))
112 >      
113 >       ntypes = getSize(atypes)
114 >      
115 >       allocate(ShapeMap%atidToShape(ntypes))
116 >    end if
117 >    
118 >    ShapeMap%currentShape = ShapeMap%currentShape + 1
119 >    current = ShapeMap%currentShape
120  
121 <    if (nAtypes == 0) then
121 >    call allocateShape(nContactFuncs, nRangeFuncs, nStrengthFuncs, &
122 >         ShapeMap%Shapes(current), stat=alloc_stat)
123 >    if (alloc_stat .ne. 0) then
124         status = -1
125         return
73 –    end if
74 –
75 –    if (.not. allocated(ShapeMap)) then
76 –       allocate(ShapeMap(nAtypes))
126      endif
127  
128 <    do i = 1, nAtypes
128 >    call getElementProperty(atypes, myAtid, "c_ident", me)
129 >    ShapeMap%atidToShape(me)                         = current
130 >    ShapeMap%Shapes(current)%atid                    = me
131 >    ShapeMap%Shapes(current)%nContactFuncs           = nContactFuncs
132 >    ShapeMap%Shapes(current)%nRangeFuncs             = nRangeFuncs
133 >    ShapeMap%Shapes(current)%nStrengthFuncs          = nStrengthFuncs
134 >    ShapeMap%Shapes(current)%ContactFuncLValue       = ContactFuncLValue
135 >    ShapeMap%Shapes(current)%ContactFuncMValue       = ContactFuncMValue
136 >    ShapeMap%Shapes(current)%ContactFunctionType     = ContactFunctionType
137 >    ShapeMap%Shapes(current)%ContactFuncCoefficient  = ContactFuncCoefficient
138 >    ShapeMap%Shapes(current)%RangeFuncLValue         = RangeFuncLValue
139 >    ShapeMap%Shapes(current)%RangeFuncMValue         = RangeFuncMValue
140 >    ShapeMap%Shapes(current)%RangeFunctionType       = RangeFunctionType
141 >    ShapeMap%Shapes(current)%RangeFuncCoefficient    = RangeFuncCoefficient
142 >    ShapeMap%Shapes(current)%StrengthFuncLValue      = StrengthFuncLValue
143 >    ShapeMap%Shapes(current)%StrengthFuncMValue      = StrengthFuncMValue
144 >    ShapeMap%Shapes(current)%StrengthFunctionType    = StrengthFunctionType
145 >    ShapeMap%Shapes(current)%StrengthFuncCoefficient = StrengthFuncCoefficient
146  
147 <       call getElementProperty(atypes, i, "is_SHAPE", thisProperty)
147 >    bigL = -1
148 >    bigM = -1
149 >    
150 >    do j = 1, ShapeMap%Shapes(current)%nContactFuncs
151 >       if (ShapeMap%Shapes(current)%ContactFuncLValue(j) .gt. bigL) then
152 >          bigL = ShapeMap%Shapes(current)%ContactFuncLValue(j)
153 >       endif
154 >       if (ShapeMap%Shapes(current)%ContactFuncMValue(j) .gt. bigM) then
155 >          bigM = ShapeMap%Shapes(current)%ContactFuncMValue(j)
156 >       endif
157 >    enddo
158 >    do j = 1, ShapeMap%Shapes(current)%nRangeFuncs
159 >       if (ShapeMap%Shapes(current)%RangeFuncLValue(j) .gt. bigL) then
160 >          bigL = ShapeMap%Shapes(current)%RangeFuncLValue(j)
161 >       endif
162 >       if (ShapeMap%Shapes(current)%RangeFuncMValue(j) .gt. bigM) then
163 >          bigM = ShapeMap%Shapes(current)%RangeFuncMValue(j)
164 >       endif
165 >    enddo
166 >    do j = 1, ShapeMap%Shapes(current)%nStrengthFuncs
167 >       if (ShapeMap%Shapes(current)%StrengthFuncLValue(j) .gt. bigL) then
168 >          bigL = ShapeMap%Shapes(current)%StrengthFuncLValue(j)
169 >       endif
170 >       if (ShapeMap%Shapes(current)%StrengthFuncMValue(j) .gt. bigM) then
171 >          bigM = ShapeMap%Shapes(current)%StrengthFuncMValue(j)
172 >       endif
173 >    enddo
174  
175 <       if (thisProperty) then
175 >    ShapeMap%Shapes(current)%bigL                    = bigL
176 >    ShapeMap%Shapes(current)%bigM                    = bigM
177  
178 <          ! do all of the shape stuff
178 >  end subroutine newShapeType
179  
180 <       endif
180 >  subroutine allocateShape(nContactFuncs, nRangeFuncs, nStrengthFuncs, &
181 >       myShape, stat)
182  
183 <       call getElementProperty(atypes, i, "is_LJ", thisProperty)
183 >    integer, intent(in) :: nContactFuncs, nRangeFuncs, nStrengthFuncs
184 >    type(Shape), intent(inout) :: myShape
185 >    integer, intent(out) :: stat
186 >    integer :: alloc_stat
187 >
188 >    if (associated(myShape%contactFuncLValue)) then
189 >       deallocate(myShape%contactFuncLValue)
190 >    endif
191 >    allocate(myShape%contactFuncLValue(nContactFuncs), stat = alloc_stat)
192 >    if (alloc_stat .ne. 0) then
193 >       stat = -1
194 >       return
195 >    endif
196 >    if (associated(myShape%contactFuncMValue)) then
197 >       deallocate(myShape%contactFuncMValue)
198 >    endif
199 >    allocate(myShape%contactFuncMValue(nContactFuncs), stat = alloc_stat)
200 >    if (alloc_stat .ne. 0) then
201 >       stat = -1
202 >       return
203 >    endif
204 >    if (associated(myShape%contactFunctionType)) then
205 >       deallocate(myShape%contactFunctionType)
206 >    endif
207 >    allocate(myShape%contactFunctionType(nContactFuncs), stat = alloc_stat)
208 >    if (alloc_stat .ne. 0) then
209 >       stat = -1
210 >       return
211 >    endif
212 >    if (associated(myShape%contactFuncCoefficient)) then
213 >       deallocate(myShape%contactFuncCoefficient)
214 >    endif
215 >    allocate(myShape%contactFuncCoefficient(nContactFuncs), stat = alloc_stat)
216 >    if (alloc_stat .ne. 0) then
217 >       stat = -1
218 >       return
219 >    endif
220  
221 <       if (thisProperty) then
222 <          ShapeMap(i)%isLJ = .true.
223 <          call getElementProperty(atypes, i, "lj_epsilon", thisDP)
224 <          ShapeMap(i)%epsilon = thisDP
225 <          call getElementProperty(atypes, i, "lj_sigma",   thisDP)
226 <          ShapeMap(i)%sigma = thisDP          
227 <       else
228 <          ShapeMap(i)%isLJ = .false.
229 <       endif
221 >    if (associated(myShape%rangeFuncLValue)) then
222 >       deallocate(myShape%rangeFuncLValue)
223 >    endif
224 >    allocate(myShape%rangeFuncLValue(nRangeFuncs), stat = alloc_stat)
225 >    if (alloc_stat .ne. 0) then
226 >       stat = -1
227 >       return
228 >    endif
229 >    if (associated(myShape%rangeFuncMValue)) then
230 >       deallocate(myShape%rangeFuncMValue)
231 >    endif
232 >    allocate(myShape%rangeFuncMValue(nRangeFuncs), stat = alloc_stat)
233 >    if (alloc_stat .ne. 0) then
234 >       stat = -1
235 >       return
236 >    endif
237 >    if (associated(myShape%rangeFunctionType)) then
238 >       deallocate(myShape%rangeFunctionType)
239 >    endif
240 >    allocate(myShape%rangeFunctionType(nRangeFuncs), stat = alloc_stat)
241 >    if (alloc_stat .ne. 0) then
242 >       stat = -1
243 >       return
244 >    endif
245 >    if (associated(myShape%rangeFuncCoefficient)) then
246 >       deallocate(myShape%rangeFuncCoefficient)
247 >    endif
248 >    allocate(myShape%rangeFuncCoefficient(nRangeFuncs), stat = alloc_stat)
249 >    if (alloc_stat .ne. 0) then
250 >       stat = -1
251 >       return
252 >    endif
253 >    
254 >    if (associated(myShape%strengthFuncLValue)) then
255 >       deallocate(myShape%strengthFuncLValue)
256 >    endif
257 >    allocate(myShape%strengthFuncLValue(nStrengthFuncs), stat = alloc_stat)
258 >    if (alloc_stat .ne. 0) then
259 >       stat = -1
260 >       return
261 >    endif
262 >    if (associated(myShape%strengthFuncMValue)) then
263 >       deallocate(myShape%strengthFuncMValue)
264 >    endif
265 >    allocate(myShape%strengthFuncMValue(nStrengthFuncs), stat = alloc_stat)
266 >    if (alloc_stat .ne. 0) then
267 >       stat = -1
268 >       return
269 >    endif
270 >    if (associated(myShape%strengthFunctionType)) then
271 >       deallocate(myShape%strengthFunctionType)
272 >    endif
273 >    allocate(myShape%strengthFunctionType(nStrengthFuncs), stat = alloc_stat)
274 >    if (alloc_stat .ne. 0) then
275 >       stat = -1
276 >       return
277 >    endif
278 >    if (associated(myShape%strengthFuncCoefficient)) then
279 >       deallocate(myShape%strengthFuncCoefficient)
280 >    endif
281 >    allocate(myShape%strengthFuncCoefficient(nStrengthFuncs), stat=alloc_stat)
282 >    if (alloc_stat .ne. 0) then
283 >       stat = -1
284 >       return
285 >    endif
286  
287 +  end subroutine allocateShape
288 +    
289 +  subroutine init_Shape_FF(status)
290 +    integer :: status
291 +    integer :: i, j, l, m, lm, function_type
292 +    real(kind=dp) :: bigSigma, myBigSigma, thisSigma, coeff, Phunc, spi
293 +    real(kind=dp) :: costheta, cpi, theta, Pi, phi, thisDP
294 +    integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, thisIP, current
295 +    logical :: thisProperty
296  
297 <    end do
297 >    Pi = 4.0d0 * datan(1.0d0)
298  
299 <    haveShapeMap = .true.
299 >    status = 0
300 >    if (ShapeMap%currentShape == 0) then
301 >       call handleError("init_Shape_FF", "No members in ShapeMap")
302 >       status = -1
303 >       return
304 >    end if
305  
306 <  end subroutine createShapeMap
306 >    bigSigma = 0.0d0
307 >    do i = 1, ShapeMap%currentShape
308  
309 +       ! Scan over theta and phi to
310 +       ! find the largest contact in any direction....
311  
312 <  
312 >       myBigSigma = 0.0d0
313 >
314 >       do iTheta = 0, nSteps
315 >          theta = (Pi/2.0d0)*(dble(iTheta)/dble(nSteps))
316 >          costheta = cos(theta)
317 >
318 >          call Associated_Legendre(costheta, ShapeMap%Shapes(i)%bigL, &
319 >               ShapeMap%Shapes(i)%bigM, lmax, plm_i, dlm_i)
320 >                  
321 >          do iPhi = 0, nSteps
322 >             phi = -Pi  + 2.0d0 * Pi * (dble(iPhi)/dble(nSteps))
323 >             cpi = cos(phi)
324 >             spi = sin(phi)
325 >            
326 >             call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(i)%bigM, &
327 >                  CHEBYSHEV_TN, tm_i, dtm_i)
328 >             call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(i)%bigM, &
329 >                  CHEBYSHEV_UN, um_i, dum_i)
330 >
331 >             thisSigma = 0.0d0
332 >
333 >             do lm = 1, ShapeMap%Shapes(i)%nContactFuncs                
334 >
335 >                l = ShapeMap%Shapes(i)%ContactFuncLValue(lm)
336 >                m = ShapeMap%Shapes(i)%ContactFuncMValue(lm)
337 >                coeff = ShapeMap%Shapes(i)%ContactFuncCoefficient(lm)
338 >                function_type = ShapeMap%Shapes(i)%ContactFunctionType(lm)
339 >
340 >                if ((function_type .eq. SH_COS).or.(m.eq.0)) then
341 >                   Phunc = coeff * tm_i(m)
342 >                else
343 >                   Phunc = coeff * spi * um_i(m-1)
344 >                endif
345 >                
346 >                thisSigma = thisSigma + plm_i(l,m)*Phunc
347 >             enddo
348 >
349 >             if (thisSigma.gt.myBigSigma) myBigSigma = thisSigma
350 >          enddo
351 >       enddo
352 >
353 >       if (myBigSigma.gt.bigSigma) bigSigma = myBigSigma
354 >    enddo
355 >
356 >    nAtypes = getSize(atypes)
357 >
358 >    if (nAtypes == 0) then
359 >       status = -1
360 >       return
361 >    end if
362 >
363 >    do i = 1, nAtypes
364 >      
365 >       call getElementProperty(atypes, i, "is_LJ", thisProperty)
366 >      
367 >       if (thisProperty) then
368 >          
369 >          ShapeMap%currentShape = ShapeMap%currentShape + 1
370 >          current = ShapeMap%currentShape
371 >
372 >          call getElementProperty(atypes, i, "c_ident",  thisIP)
373 >          ShapeMap%atidToShape(thisIP) = current
374 >          ShapeMap%Shapes(current)%atid = thisIP
375 >
376 >          ShapeMap%Shapes(current)%isLJ = .true.
377 >
378 >          call getElementProperty(atypes, i, "lj_epsilon", thisDP)
379 >          ShapeMap%Shapes(current)%epsilon = thisDP
380 >
381 >          call getElementProperty(atypes, i, "lj_sigma",   thisDP)
382 >          ShapeMap%Shapes(current)%sigma = thisDP          
383 >          if (thisDP .gt. bigSigma) bigSigma = thisDP
384 >          
385 >       endif
386 >      
387 >    end do
388 >
389 >    haveShapeMap = .true.
390 >    
391 >  end subroutine init_Shape_FF
392 >    
393    subroutine do_shape_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
394         pot, A, f, t, do_pot)
395      
# Line 121 | Line 404 | contains  
404      logical, intent(in) :: do_pot
405  
406      real (kind=dp) :: r3, r5, rt2, rt3, rt5, rt6, rt11, rt12, rt126
407 <    integer :: me1, me2, l, m, lm, id1, id2, localError, function_type
407 >    integer :: atid1, atid2, st1, st2
408 >    integer :: l, m, lm, id1, id2, localError, function_type
409      real (kind=dp) :: sigma_i, s_i, eps_i, sigma_j, s_j, eps_j
410      real (kind=dp) :: coeff
411  
# Line 201 | Line 485 | contains  
485      real (kind=dp) :: fxradial, fyradial, fzradial
486  
487      if (.not.haveShapeMap) then
488 <       localError = 0
489 <       call createShapeMap(localError)
206 <       if ( localError .ne. 0 ) then
207 <          call handleError("calc_shape", "ShapeMap creation failed!")
208 <          return
209 <       end if
488 >       call handleError("calc_shape", "NO SHAPEMAP!!!!")
489 >       return      
490      endif
491      
492      !! We assume that the rotation matrices have already been calculated
# Line 223 | Line 503 | contains  
503      drdyj = d(2) / rij
504      drdzj = d(3) / rij
505      
506 +    ! find the atom type id (atid) for each atom:
507   #ifdef IS_MPI
508 <    me1 = atid_Row(atom1)
509 <    me2 = atid_Col(atom2)
508 >    atid1 = atid_Row(atom1)
509 >    atid2 = atid_Col(atom2)
510   #else
511 <    me1 = atid(atom1)
512 <    me2 = atid(atom2)
511 >    atid1 = atid(atom1)
512 >    atid2 = atid(atom2)
513   #endif
514  
515 <    if (ShapeMap(me1)%isLJ) then
516 <       sigma_i = ShapeMap(me1)%sigma
517 <       s_i = ShapeMap(me1)%sigma
518 <       eps_i = ShapeMap(me1)%epsilon
515 >    ! use the atid to find the shape type (st) for each atom:
516 >
517 >    st1 = ShapeMap%atidToShape(atid1)
518 >    st2 = ShapeMap%atidToShape(atid2)
519 >    
520 >    if (ShapeMap%Shapes(st1)%isLJ) then
521 >       sigma_i = ShapeMap%Shapes(st1)%sigma
522 >       s_i = ShapeMap%Shapes(st1)%sigma
523 >       eps_i = ShapeMap%Shapes(st1)%epsilon
524         dsigmaidx = 0.0d0
525         dsigmaidy = 0.0d0
526         dsigmaidz = 0.0d0
# Line 304 | Line 590 | contains  
590         dspiduy = xi * yi * zi / proji3
591         dspiduz = xi * (1.0d0 / proji - yi2 / proji3) + (xi * yi2 / proji3)
592  
593 <       call Associated_Legendre(cti, ShapeMap(me1)%bigL, &
594 <            ShapeMap(me1)%bigM, lmax, plm_i, dlm_i)
593 >       call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigL, &
594 >            ShapeMap%Shapes(st1)%bigM, lmax, plm_i, dlm_i)
595  
596 <       call Orthogonal_Polynomial(cpi, ShapeMap(me1)%bigM, CHEBYSHEV_TN, &
597 <            tm_i, dtm_i)
598 <       call Orthogonal_Polynomial(cpi, ShapeMap(me1)%bigM, CHEBYSHEV_UN, &
599 <            um_i, dum_i)
596 >       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
597 >            CHEBYSHEV_TN, tm_i, dtm_i)
598 >       call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
599 >            CHEBYSHEV_UN, um_i, dum_i)
600        
601         sigma_i = 0.0d0
602         s_i = 0.0d0
# Line 334 | Line 620 | contains  
620         depsiduy = 0.0d0
621         depsiduz = 0.0d0
622  
623 <       do lm = 1, ShapeMap(me1)%nContactFuncs
624 <          l = ShapeMap(me1)%ContactFuncLValue(lm)
625 <          m = ShapeMap(me1)%ContactFuncMValue(lm)
626 <          coeff = ShapeMap(me1)%ContactFuncCoefficient(lm)
627 <          function_type = ShapeMap(me1)%ContactFunctionType(lm)
623 >       do lm = 1, ShapeMap%Shapes(st1)%nContactFuncs
624 >          l = ShapeMap%Shapes(st1)%ContactFuncLValue(lm)
625 >          m = ShapeMap%Shapes(st1)%ContactFuncMValue(lm)
626 >          coeff = ShapeMap%Shapes(st1)%ContactFuncCoefficient(lm)
627 >          function_type = ShapeMap%Shapes(st1)%ContactFunctionType(lm)
628  
629            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
630               Phunc = coeff * tm_i(m)
# Line 376 | Line 662 | contains  
662  
663         end do
664  
665 <       do lm = 1, ShapeMap(me1)%nRangeFuncs
666 <          l = ShapeMap(me1)%RangeFuncLValue(lm)
667 <          m = ShapeMap(me1)%RangeFuncMValue(lm)
668 <          coeff = ShapeMap(me1)%RangeFuncCoefficient(lm)
669 <          function_type = ShapeMap(me1)%RangeFunctionType(lm)
665 >       do lm = 1, ShapeMap%Shapes(st1)%nRangeFuncs
666 >          l = ShapeMap%Shapes(st1)%RangeFuncLValue(lm)
667 >          m = ShapeMap%Shapes(st1)%RangeFuncMValue(lm)
668 >          coeff = ShapeMap%Shapes(st1)%RangeFuncCoefficient(lm)
669 >          function_type = ShapeMap%Shapes(st1)%RangeFunctionType(lm)
670            
671            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
672               Phunc = coeff * tm_i(m)
# Line 418 | Line 704 | contains  
704  
705         end do
706                
707 <       do lm = 1, ShapeMap(me1)%nStrengthFuncs
708 <          l = ShapeMap(me1)%StrengthFuncLValue(lm)
709 <          m = ShapeMap(me1)%StrengthFuncMValue(lm)
710 <          coeff = ShapeMap(me1)%StrengthFuncCoefficient(lm)
711 <          function_type = ShapeMap(me1)%StrengthFunctionType(lm)
707 >       do lm = 1, ShapeMap%Shapes(st1)%nStrengthFuncs
708 >          l = ShapeMap%Shapes(st1)%StrengthFuncLValue(lm)
709 >          m = ShapeMap%Shapes(st1)%StrengthFuncMValue(lm)
710 >          coeff = ShapeMap%Shapes(st1)%StrengthFuncCoefficient(lm)
711 >          function_type = ShapeMap%Shapes(st1)%StrengthFunctionType(lm)
712            
713            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
714               Phunc = coeff * tm_i(m)
# Line 464 | Line 750 | contains  
750        
751         ! now do j:
752  
753 <    if (ShapeMap(me2)%isLJ) then
754 <       sigma_j = ShapeMap(me2)%sigma
755 <       s_j = ShapeMap(me2)%sigma
756 <       eps_j = ShapeMap(me2)%epsilon
753 >    if (ShapeMap%Shapes(st2)%isLJ) then
754 >       sigma_j = ShapeMap%Shapes(st2)%sigma
755 >       s_j = ShapeMap%Shapes(st2)%sigma
756 >       eps_j = ShapeMap%Shapes(st2)%epsilon
757         dsigmajdx = 0.0d0
758         dsigmajdy = 0.0d0
759         dsigmajdz = 0.0d0
# Line 537 | Line 823 | contains  
823         dspjduy = xj * yj * zj / projj3
824         dspjduz = xj * (1.0d0 / projj - yi2 / projj3) + (xj * yj2 / projj3)
825        
826 <       call Associated_Legendre(ctj, ShapeMap(me2)%bigL, &
827 <            ShapeMap(me2)%bigM, lmax, plm_j, dlm_j)
826 >       call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigL, &
827 >            ShapeMap%Shapes(st2)%bigM, lmax, plm_j, dlm_j)
828        
829 <       call Orthogonal_Polynomial(cpj, ShapeMap(me2)%bigM, CHEBYSHEV_TN, &
830 <            tm_j, dtm_j)
831 <       call Orthogonal_Polynomial(cpj, ShapeMap(me2)%bigM, CHEBYSHEV_UN, &
832 <            um_j, dum_j)
829 >       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
830 >            CHEBYSHEV_TN, tm_j, dtm_j)
831 >       call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
832 >            CHEBYSHEV_UN, um_j, dum_j)
833        
834         sigma_j = 0.0d0
835         s_j = 0.0d0
# Line 567 | Line 853 | contains  
853         depsjduy = 0.0d0
854         depsjduz = 0.0d0
855  
856 <       do lm = 1, ShapeMap(me2)%nContactFuncs
857 <          l = ShapeMap(me2)%ContactFuncLValue(lm)
858 <          m = ShapeMap(me2)%ContactFuncMValue(lm)
859 <          coeff = ShapeMap(me2)%ContactFuncCoefficient(lm)
860 <          function_type = ShapeMap(me2)%ContactFunctionType(lm)
856 >       do lm = 1, ShapeMap%Shapes(st2)%nContactFuncs
857 >          l = ShapeMap%Shapes(st2)%ContactFuncLValue(lm)
858 >          m = ShapeMap%Shapes(st2)%ContactFuncMValue(lm)
859 >          coeff = ShapeMap%Shapes(st2)%ContactFuncCoefficient(lm)
860 >          function_type = ShapeMap%Shapes(st2)%ContactFunctionType(lm)
861  
862            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
863               Phunc = coeff * tm_j(m)
# Line 609 | Line 895 | contains  
895  
896         end do
897  
898 <       do lm = 1, ShapeMap(me2)%nRangeFuncs
899 <          l = ShapeMap(me2)%RangeFuncLValue(lm)
900 <          m = ShapeMap(me2)%RangeFuncMValue(lm)
901 <          coeff = ShapeMap(me2)%RangeFuncCoefficient(lm)
902 <          function_type = ShapeMap(me2)%RangeFunctionType(lm)
898 >       do lm = 1, ShapeMap%Shapes(st2)%nRangeFuncs
899 >          l = ShapeMap%Shapes(st2)%RangeFuncLValue(lm)
900 >          m = ShapeMap%Shapes(st2)%RangeFuncMValue(lm)
901 >          coeff = ShapeMap%Shapes(st2)%RangeFuncCoefficient(lm)
902 >          function_type = ShapeMap%Shapes(st2)%RangeFunctionType(lm)
903  
904            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
905               Phunc = coeff * tm_j(m)
# Line 651 | Line 937 | contains  
937  
938         end do
939  
940 <       do lm = 1, ShapeMap(me2)%nStrengthFuncs
941 <          l = ShapeMap(me2)%StrengthFuncLValue(lm)
942 <          m = ShapeMap(me2)%StrengthFuncMValue(lm)
943 <          coeff = ShapeMap(me2)%StrengthFuncCoefficient(lm)
944 <          function_type = ShapeMap(me2)%StrengthFunctionType(lm)
940 >       do lm = 1, ShapeMap%Shapes(st2)%nStrengthFuncs
941 >          l = ShapeMap%Shapes(st2)%StrengthFuncLValue(lm)
942 >          m = ShapeMap%Shapes(st2)%StrengthFuncMValue(lm)
943 >          coeff = ShapeMap%Shapes(st2)%StrengthFuncCoefficient(lm)
944 >          function_type = ShapeMap%Shapes(st2)%StrengthFunctionType(lm)
945  
946            if ((function_type .eq. SH_COS).or.(m.eq.0)) then
947               Phunc = coeff * tm_j(m)

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines