ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/OOPSE-4/src/UseTheForce/DarkSide/shapes.F90

Comparing trunk/OOPSE-4/src/UseTheForce/DarkSide/shapes.F90 (file contents):
Revision 1633 by gezelter, Fri Oct 22 20:22:48 2004 UTC vs.
Revision 2357 by chuckv, Wed Oct 12 20:18:17 2005 UTC

#	Line 1 | Line 1
1	+	!!
2	+	!! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	+	!!
4	+	!! The University of Notre Dame grants you ("Licensee") a
5	+	!! non-exclusive, royalty free, license to use, modify and
6	+	!! redistribute this software in source and binary code form, provided
7	+	!! that the following conditions are met:
8	+	!!
9	+	!! 1. Acknowledgement of the program authors must be made in any
10	+	!!    publication of scientific results based in part on use of the
11	+	!!    program.  An acceptable form of acknowledgement is citation of
12	+	!!    the article in which the program was described (Matthew
13	+	!!    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	+	!!    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	+	!!    Parallel Simulation Engine for Molecular Dynamics,"
16	+	!!    J. Comput. Chem. 26, pp. 252-271 (2005))
17	+	!!
18	+	!! 2. Redistributions of source code must retain the above copyright
19	+	!!    notice, this list of conditions and the following disclaimer.
20	+	!!
21	+	!! 3. Redistributions in binary form must reproduce the above copyright
22	+	!!    notice, this list of conditions and the following disclaimer in the
23	+	!!    documentation and/or other materials provided with the
24	+	!!    distribution.
25	+	!!
26	+	!! This software is provided "AS IS," without a warranty of any
27	+	!! kind. All express or implied conditions, representations and
28	+	!! warranties, including any implied warranty of merchantability,
29	+	!! fitness for a particular purpose or non-infringement, are hereby
30	+	!! excluded.  The University of Notre Dame and its licensors shall not
31	+	!! be liable for any damages suffered by licensee as a result of
32	+	!! using, modifying or distributing the software or its
33	+	!! derivatives. In no event will the University of Notre Dame or its
34	+	!! licensors be liable for any lost revenue, profit or data, or for
35	+	!! direct, indirect, special, consequential, incidental or punitive
36	+	!! damages, however caused and regardless of the theory of liability,
37	+	!! arising out of the use of or inability to use software, even if the
38	+	!! University of Notre Dame has been advised of the possibility of
39	+	!! such damages.
40	+	!!
41	+
42	+
43		module shapes
44
45		use force_globals
#	Line 13 | Line 55 | module shapes
55		implicit none
56
57		PRIVATE
58	<
58	>	#define __FORTRAN90
59	>	#include "UseTheForce/DarkSide/fInteractionMap.h"
60		INTEGER, PARAMETER:: CHEBYSHEV_TN = 1
61		INTEGER, PARAMETER:: CHEBYSHEV_UN = 2
62		INTEGER, PARAMETER:: LAGUERRE     = 3
#	Line 25 | Line 68 | module shapes
68
69		public :: do_shape_pair
70		public :: newShapeType
71	+	public :: complete_Shape_FF
72	+	public :: destroyShapeTypes
73	+	public :: getShapeCut
74
29	–
75		type, private :: Shape
76		integer :: atid
77		integer :: nContactFuncs
#	Line 50 | Line 95 | module shapes
95		real ( kind = dp )  :: epsilon
96		real ( kind = dp )  :: sigma
97		end type Shape
98	<
98	>
99		type, private :: ShapeList
100		integer :: n_shapes = 0
101		integer :: currentShape = 0
102	<	type (Shape), pointer :: Shapes(:)      => null()
103	<	integer, pointer      :: atidToShape(:) => null()
102	>	type(Shape), pointer :: Shapes(:)      => null()
103	>	integer, pointer     :: atidToShape(:) => null()
104		end type ShapeList
105
106		type(ShapeList), save :: ShapeMap
107	<
107	>
108		integer :: lmax
109	<	real (kind=dp), allocatable, dimension(:,:) :: plm_i, dlm_i, plm_j, dlm_j
65	<	real (kind=dp), allocatable, dimension(:) :: tm_i, dtm_i, um_i, dum_i
66	<	real (kind=dp), allocatable, dimension(:) :: tm_j, dtm_j, um_j, dum_j
67	<
109	>
110		contains
111	<
111	>
112		subroutine newShapeType(nContactFuncs, ContactFuncLValue, &
113		ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
114		nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
115		RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
116		StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
117	<	myAtid, status)
118	<
117	>	c_ident, status)
118	>
119		integer :: nContactFuncs
120		integer :: nRangeFuncs
121		integer :: nStrengthFuncs
122		integer :: shape_ident
123		integer :: status
124	<	integer :: myAtid
124	>	integer :: c_ident
125	>	integer :: myATID
126		integer :: bigL
127		integer :: bigM
128		integer :: j, me, nShapeTypes, nLJTypes, ntypes, current, alloc_stat
#	Line 104 | Line 147 | contains
147
148		call getMatchingElementList(atypes, "is_Shape", .true., &
149		nShapeTypes, MatchList)
150	<
150	>
151		call getMatchingElementList(atypes, "is_LennardJones", .true., &
152		nLJTypes, MatchList)
153	<
153	>
154		ShapeMap%n_shapes = nShapeTypes + nLJTypes
155	<
155	>
156		allocate(ShapeMap%Shapes(nShapeTypes + nLJTypes))
157	<
157	>
158		ntypes = getSize(atypes)
159	<
159	>
160		allocate(ShapeMap%atidToShape(ntypes))
161		end if
162	<
162	>
163		ShapeMap%currentShape = ShapeMap%currentShape + 1
164		current = ShapeMap%currentShape
165
#	Line 127 | Line 170 | contains
170		return
171		endif
172
173	<	call getElementProperty(atypes, myAtid, "c_ident", me)
174	<	ShapeMap%atidToShape(me)                         = current
175	<	ShapeMap%Shapes(current)%atid                    = me
173	>	myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
174	>
175	>	ShapeMap%atidToShape(myATID)                     = current
176	>	ShapeMap%Shapes(current)%atid                    = myATID
177		ShapeMap%Shapes(current)%nContactFuncs           = nContactFuncs
178		ShapeMap%Shapes(current)%nRangeFuncs             = nRangeFuncs
179		ShapeMap%Shapes(current)%nStrengthFuncs          = nStrengthFuncs
#	Line 148 | Line 192 | contains
192
193		bigL = -1
194		bigM = -1
195	<
195	>
196		do j = 1, ShapeMap%Shapes(current)%nContactFuncs
197		if (ShapeMap%Shapes(current)%ContactFuncLValue(j) .gt. bigL) then
198		bigL = ShapeMap%Shapes(current)%ContactFuncLValue(j)
#	Line 186 | Line 230 | contains
230		type(Shape), intent(inout) :: myShape
231		integer, intent(out) :: stat
232		integer :: alloc_stat
233	<
233	>
234	>	stat = 0
235		if (associated(myShape%contactFuncLValue)) then
236		deallocate(myShape%contactFuncLValue)
237		endif
#	Line 252 | Line 297 | contains
297		stat = -1
298		return
299		endif
300	<
300	>
301		if (associated(myShape%strengthFuncLValue)) then
302		deallocate(myShape%strengthFuncLValue)
303		endif
#	Line 286 | Line 331 | contains
331		return
332		endif
333
334	+	return
335	+
336		end subroutine allocateShape
337	<
338	<	subroutine init_Shape_FF(status)
337	>
338	>	subroutine complete_Shape_FF(status)
339		integer :: status
340		integer :: i, j, l, m, lm, function_type
341	<	real(kind=dp) :: bigSigma, myBigSigma, thisSigma, coeff, Phunc, spi
342	<	real(kind=dp) :: costheta, cpi, theta, Pi, phi, thisDP, sigma
296	<	integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, thisIP, current
341	>	real(kind=dp) :: thisDP, sigma
342	>	integer :: alloc_stat, iTheta, iPhi, nSteps, nAtypes, myATID, current
343		logical :: thisProperty
344
299	–	Pi = 4.0d0 * datan(1.0d0)
300	–
345		status = 0
346		if (ShapeMap%currentShape == 0) then
347		call handleError("init_Shape_FF", "No members in ShapeMap")
348		status = -1
349		return
350		end if
307	–
308	–	bigSigma = 0.0d0
309	–	do i = 1, ShapeMap%currentShape
310	–
311	–	! Scan over theta and phi to
312	–	! find the largest contact in any direction....
313	–
314	–	myBigSigma = 0.0d0
315	–
316	–	do iTheta = 0, nSteps
317	–	theta = (Pi/2.0d0)*(dble(iTheta)/dble(nSteps))
318	–	costheta = cos(theta)
319	–
320	–	call Associated_Legendre(costheta, ShapeMap%Shapes(i)%bigL, &
321	–	ShapeMap%Shapes(i)%bigM, lmax, plm_i, dlm_i)
322	–
323	–	do iPhi = 0, nSteps
324	–	phi = -Pi  + 2.0d0 * Pi * (dble(iPhi)/dble(nSteps))
325	–	cpi = cos(phi)
326	–	spi = sin(phi)
327	–
328	–	call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(i)%bigM, &
329	–	CHEBYSHEV_TN, tm_i, dtm_i)
330	–	call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(i)%bigM, &
331	–	CHEBYSHEV_UN, um_i, dum_i)
351
333	–	thisSigma = 0.0d0
334	–
335	–	do lm = 1, ShapeMap%Shapes(i)%nContactFuncs
336	–
337	–	l = ShapeMap%Shapes(i)%ContactFuncLValue(lm)
338	–	m = ShapeMap%Shapes(i)%ContactFuncMValue(lm)
339	–	coeff = ShapeMap%Shapes(i)%ContactFuncCoefficient(lm)
340	–	function_type = ShapeMap%Shapes(i)%ContactFunctionType(lm)
341	–
342	–	if ((function_type .eq. SH_COS).or.(m.eq.0)) then
343	–	Phunc = coeff * tm_i(m)
344	–	else
345	–	Phunc = coeff * spi * um_i(m-1)
346	–	endif
347	–
348	–	thisSigma = thisSigma + plm_i(l,m)*Phunc
349	–	enddo
350	–
351	–	if (thisSigma.gt.myBigSigma) myBigSigma = thisSigma
352	–	enddo
353	–	enddo
354	–
355	–	if (myBigSigma.gt.bigSigma) bigSigma = myBigSigma
356	–	enddo
357	–
352		nAtypes = getSize(atypes)
353
354		if (nAtypes == 0) then
355		status = -1
356		return
357	<	end if
357	>	end if
358
359	+	! atypes comes from c side
360		do i = 1, nAtypes
361	<
362	<	call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
363	<
361	>
362	>	myATID = getFirstMatchingElement(atypes, 'c_ident', i)
363	>	call getElementProperty(atypes, myATID, "is_LennardJones", thisProperty)
364	>
365		if (thisProperty) then
370	–
366		ShapeMap%currentShape = ShapeMap%currentShape + 1
367		current = ShapeMap%currentShape
368
369	<	call getElementProperty(atypes, i, "c_ident",  thisIP)
370	<	ShapeMap%atidToShape(thisIP) = current
376	<	ShapeMap%Shapes(current)%atid = thisIP
369	>	ShapeMap%atidToShape(myATID) = current
370	>	ShapeMap%Shapes(current)%atid = myATID
371
372		ShapeMap%Shapes(current)%isLJ = .true.
373
374	<	ShapeMap%Shapes(current)%epsilon = getEpsilon(thisIP)
375	<	sigma = getSigma(thisIP)
376	<	ShapeMap%Shapes(current)%sigma = sigma
383	<	if (sigma .gt. bigSigma) bigSigma = thisDP
384	<
374	>	ShapeMap%Shapes(current)%epsilon = getEpsilon(myATID)
375	>	ShapeMap%Shapes(current)%sigma = getSigma(myATID)
376	>
377		endif
378	<
378	>
379		end do
380
381		haveShapeMap = .true.
382	<
383	<	end subroutine init_Shape_FF
384	<
382	>
383	>	!    do i = 1, ShapeMap%n_shapes
384	>	!       write(*,*) 'i = ', i, ' isLJ = ', ShapeMap%Shapes(i)%isLJ
385	>	!    end do
386	>
387	>	end subroutine complete_Shape_FF
388	>
389	>	function getShapeCut(atomID) result(cutValue)
390	>	integer, intent(in) :: atomID
391	>	real(kind=dp) :: cutValue, whoopdedoo
392	>
393	>	!! this is just a placeholder for a cutoff value, hopefully we'll
394	>	!! develop a method to calculate a sensible value
395	>	whoopdedoo = 9.0_dp
396	>
397	>	cutValue = whoopdedoo
398	>
399	>	end function getShapeCut
400	>
401		subroutine do_shape_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, &
402		pot, A, f, t, do_pot)
403	<
403	>
404	>	INTEGER, PARAMETER:: LMAX         = 64
405	>	INTEGER, PARAMETER:: MMAX         = 64
406	>
407		integer, intent(in) :: atom1, atom2
408		real (kind=dp), intent(inout) :: rij, r2
409		real (kind=dp), dimension(3), intent(in) :: d
410		real (kind=dp), dimension(3), intent(inout) :: fpair
411	<	real (kind=dp) :: pot, vpair, sw
411	>	real (kind=dp) :: pot, vpair, sw, dswdr
412		real (kind=dp), dimension(9,nLocal) :: A
413		real (kind=dp), dimension(3,nLocal) :: f
414		real (kind=dp), dimension(3,nLocal) :: t
#	Line 408 | Line 419 | contains
419		integer :: l, m, lm, id1, id2, localError, function_type
420		real (kind=dp) :: sigma_i, s_i, eps_i, sigma_j, s_j, eps_j
421		real (kind=dp) :: coeff
422	+	real (kind=dp) :: pot_temp
423
424		real (kind=dp) :: dsigmaidx, dsigmaidy, dsigmaidz
425		real (kind=dp) :: dsigmaidux, dsigmaiduy, dsigmaiduz
#	Line 426 | Line 438 | contains
438
439		real (kind=dp) :: xi, yi, zi, xj, yj, zj, xi2, yi2, zi2, xj2, yj2, zj2
440
441	+	real (kind=dp) :: sti2, stj2
442	+
443		real (kind=dp) :: proji, proji3, projj, projj3
444		real (kind=dp) :: cti, ctj, cpi, cpj, spi, spj
445		real (kind=dp) :: Phunc, sigma, s, eps, rtdenom, rt
#	Line 457 | Line 471 | contains
471		real (kind=dp) :: dsduxi, dsduyi, dsduzi
472		real (kind=dp) :: dsdxj, dsdyj, dsdzj
473		real (kind=dp) :: dsduxj, dsduyj, dsduzj
474	<
474	>
475		real (kind=dp) :: depsdxi, depsdyi, depsdzi
476		real (kind=dp) :: depsduxi, depsduyi, depsduzi
477		real (kind=dp) :: depsdxj, depsdyj, depsdzj
#	Line 484 | Line 498 | contains
498		real (kind=dp) :: fxji, fyji, fzji, fxjj, fyjj, fzjj
499		real (kind=dp) :: fxradial, fyradial, fzradial
500
501	+	real (kind=dp) :: xihat, yihat, zihat, xjhat, yjhat, zjhat
502	+
503	+	real (kind=dp) :: plm_i(0:LMAX,0:MMAX), dlm_i(0:LMAX,0:MMAX)
504	+	real (kind=dp) :: plm_j(0:LMAX,0:MMAX), dlm_j(0:LMAX,0:MMAX)
505	+	real (kind=dp) :: tm_i(0:MMAX), dtm_i(0:MMAX), um_i(0:MMAX), dum_i(0:MMAX)
506	+	real (kind=dp) :: tm_j(0:MMAX), dtm_j(0:MMAX), um_j(0:MMAX), dum_j(0:MMAX)
507	+
508		if (.not.haveShapeMap) then
509		call handleError("calc_shape", "NO SHAPEMAP!!!!")
510		return
511		endif
512	<
512	>
513		!! We assume that the rotation matrices have already been calculated
514		!! and placed in the A array.
494	–
515		r3 = r2*rij
516		r5 = r3*r2
517	<
517	>
518		drdxi = -d(1) / rij
519		drdyi = -d(2) / rij
520		drdzi = -d(3) / rij
521	+	drduxi = 0.0d0
522	+	drduyi = 0.0d0
523	+	drduzi = 0.0d0
524
525		drdxj = d(1) / rij
526		drdyj = d(2) / rij
527		drdzj = d(3) / rij
528	<
528	>	drduxj = 0.0d0
529	>	drduyj = 0.0d0
530	>	drduzj = 0.0d0
531	>
532		! find the atom type id (atid) for each atom:
533		#ifdef IS_MPI
534		atid1 = atid_Row(atom1)
#	Line 513 | Line 539 | contains
539		#endif
540
541		! use the atid to find the shape type (st) for each atom:
516	–
542		st1 = ShapeMap%atidToShape(atid1)
543		st2 = ShapeMap%atidToShape(atid2)
544
545	+	!    write(*,*) atom1, atom2, atid1, atid2, st1, st2, ShapeMap%Shapes(st1)%isLJ, ShapeMap%Shapes(st2)%isLJ
546	+
547		if (ShapeMap%Shapes(st1)%isLJ) then
548	+
549		sigma_i = ShapeMap%Shapes(st1)%sigma
550		s_i = ShapeMap%Shapes(st1)%sigma
551		eps_i = ShapeMap%Shapes(st1)%epsilon
#	Line 544 | Line 572 | contains
572		#ifdef IS_MPI
573		! rotate the inter-particle separation into the two different
574		! body-fixed coordinate systems:
575	<
575	>
576		xi = A_row(1,atom1)*d(1) + A_row(2,atom1)*d(2) + A_row(3,atom1)*d(3)
577		yi = A_row(4,atom1)*d(1) + A_row(5,atom1)*d(2) + A_row(6,atom1)*d(3)
578		zi = A_row(7,atom1)*d(1) + A_row(8,atom1)*d(2) + A_row(9,atom1)*d(3)
579	<
579	>
580		#else
581		! rotate the inter-particle separation into the two different
582		! body-fixed coordinate systems:
583	<
583	>
584		xi = a(1,atom1)*d(1) + a(2,atom1)*d(2) + a(3,atom1)*d(3)
585		yi = a(4,atom1)*d(1) + a(5,atom1)*d(2) + a(6,atom1)*d(3)
586		zi = a(7,atom1)*d(1) + a(8,atom1)*d(2) + a(9,atom1)*d(3)
587	<
587	>
588		#endif
589	<
589	>	xihat = xi / rij
590	>	yihat = yi / rij
591	>	zihat = zi / rij
592		xi2 = xi*xi
593		yi2 = yi*yi
594	<	zi2 = zi*zi
565	<
566	<	proji = sqrt(xi2 + yi2)
567	<	proji3 = proji*proji*proji
568	<
594	>	zi2 = zi*zi
595		cti = zi / rij
596	+
597	+	if (cti .gt. 1.0_dp) cti = 1.0_dp
598	+	if (cti .lt. -1.0_dp) cti = -1.0_dp
599	+
600		dctidx = - zi * xi / r3
601		dctidy = - zi * yi / r3
602		dctidz = 1.0d0 / rij - zi2 / r3
603	<	dctidux =  yi / rij
604	<	dctiduy = -xi / rij
605	<	dctiduz = 0.0d0
606	<
603	>	dctidux = yi / rij ! - (zi * xi2) / r3
604	>	dctiduy = -xi / rij !- (zi * yi2) / r3
605	>	dctiduz = 0.0d0 !zi / rij - (zi2 * zi) / r3
606	>
607	>	! this is an attempt to try to truncate the singularity when
608	>	! sin(theta) is near 0.0:
609	>
610	>	sti2 = 1.0_dp - cti*cti
611	>	if (dabs(sti2) .lt. 1.0d-12) then
612	>	proji = sqrt(rij * 1.0d-12)
613	>	dcpidx = 1.0d0 / proji
614	>	dcpidy = 0.0d0
615	>	dcpidux = xi / proji
616	>	dcpiduy = 0.0d0
617	>	dspidx = 0.0d0
618	>	dspidy = 1.0d0 / proji
619	>	dspidux = 0.0d0
620	>	dspiduy = yi / proji
621	>	else
622	>	proji = sqrt(xi2 + yi2)
623	>	proji3 = proji*proji*proji
624	>	dcpidx = 1.0d0 / proji - xi2 / proji3
625	>	dcpidy = - xi * yi / proji3
626	>	dcpidux = xi / proji - (xi2 * xi) / proji3
627	>	dcpiduy = - (xi * yi2) / proji3
628	>	dspidx = - xi * yi / proji3
629	>	dspidy = 1.0d0 / proji - yi2 / proji3
630	>	dspidux = - (yi * xi2) / proji3
631	>	dspiduy = yi / proji - (yi2 * yi) / proji3
632	>	endif
633	>
634		cpi = xi / proji
578	–	dcpidx = 1.0d0 / proji - xi2 / proji3
579	–	dcpidy = - xi * yi / proji3
635		dcpidz = 0.0d0
636	<	dcpidux = xi * yi * zi / proji3
637	<	dcpiduy = -zi * (1.0d0 / proji - xi2 / proji3)
583	<	dcpiduz = -yi * (1.0d0 / proji - xi2 / proji3)  - (xi2 * yi / proji3)
584	<
636	>	dcpiduz = 0.0d0
637	>
638		spi = yi / proji
586	–	dspidx = - xi * yi / proji3
587	–	dspidy = 1.0d0 / proji - yi2 / proji3
639		dspidz = 0.0d0
640	<	dspidux = -zi * (1.0d0 / proji - yi2 / proji3)
590	<	dspiduy = xi * yi * zi / proji3
591	<	dspiduz = xi * (1.0d0 / proji - yi2 / proji3) + (xi * yi2 / proji3)
640	>	dspiduz = 0.0d0
641
642	<	call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigL, &
643	<	ShapeMap%Shapes(st1)%bigM, lmax, plm_i, dlm_i)
642	>	call Associated_Legendre(cti, ShapeMap%Shapes(st1)%bigM, &
643	>	ShapeMap%Shapes(st1)%bigL, LMAX, &
644	>	plm_i, dlm_i)
645
646	<	call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
646	>	call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
647		CHEBYSHEV_TN, tm_i, dtm_i)
648	<	call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, &
648	>	call Orthogonal_Polynomial(cpi, ShapeMap%Shapes(st1)%bigM, MMAX, &
649		CHEBYSHEV_UN, um_i, dum_i)
650	<
650	>
651		sigma_i = 0.0d0
652		s_i = 0.0d0
653		eps_i = 0.0d0
#	Line 631 | Line 681 | contains
681		dPhuncdX = coeff * dtm_i(m) * dcpidx
682		dPhuncdY = coeff * dtm_i(m) * dcpidy
683		dPhuncdZ = coeff * dtm_i(m) * dcpidz
684	<	dPhuncdUz = coeff * dtm_i(m) * dcpidux
684	>	dPhuncdUx = coeff * dtm_i(m) * dcpidux
685		dPhuncdUy = coeff * dtm_i(m) * dcpiduy
686		dPhuncdUz = coeff * dtm_i(m) * dcpiduz
687		else
#	Line 644 | Line 694 | contains
694		dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
695		endif
696
697	<	sigma_i = sigma_i + plm_i(l,m)*Phunc
698	<
699	<	dsigmaidx = dsigmaidx + plm_i(l,m)*dPhuncdX + &
700	<	Phunc * dlm_i(l,m) * dctidx
701	<	dsigmaidy = dsigmaidy + plm_i(l,m)*dPhuncdY + &
702	<	Phunc * dlm_i(l,m) * dctidy
703	<	dsigmaidz = dsigmaidz + plm_i(l,m)*dPhuncdZ + &
704	<	Phunc * dlm_i(l,m) * dctidz
705	<
706	<	dsigmaidux = dsigmaidux + plm_i(l,m)* dPhuncdUx + &
707	<	Phunc * dlm_i(l,m) * dctidux
708	<	dsigmaiduy = dsigmaiduy + plm_i(l,m)* dPhuncdUy + &
709	<	Phunc * dlm_i(l,m) * dctiduy
710	<	dsigmaiduz = dsigmaiduz + plm_i(l,m)* dPhuncdUz + &
711	<	Phunc * dlm_i(l,m) * dctiduz
712	<
697	>	sigma_i = sigma_i + plm_i(m,l)*Phunc
698	>	!!$          write(*,*) 'dsigmaidux = ', dsigmaidux
699	>	!!$          write(*,*) 'Phunc = ', Phunc
700	>	dsigmaidx = dsigmaidx + plm_i(m,l)*dPhuncdX + &
701	>	Phunc * dlm_i(m,l) * dctidx
702	>	dsigmaidy = dsigmaidy + plm_i(m,l)*dPhuncdY + &
703	>	Phunc * dlm_i(m,l) * dctidy
704	>	dsigmaidz = dsigmaidz + plm_i(m,l)*dPhuncdZ + &
705	>	Phunc * dlm_i(m,l) * dctidz
706	>	dsigmaidux = dsigmaidux + plm_i(m,l)* dPhuncdUx + &
707	>	Phunc * dlm_i(m,l) * dctidux
708	>	dsigmaiduy = dsigmaiduy + plm_i(m,l)* dPhuncdUy + &
709	>	Phunc * dlm_i(m,l) * dctiduy
710	>	dsigmaiduz = dsigmaiduz + plm_i(m,l)* dPhuncdUz + &
711	>	Phunc * dlm_i(m,l) * dctiduz
712	>	!!$          write(*,*) 'dsigmaidux = ', dsigmaidux, '; dPhuncdUx = ', dPhuncdUx, &
713	>	!!$                     '; dctidux = ', dctidux, '; plm_i(m,l) = ', plm_i(m,l), &
714	>	!!$                     '; dlm_i(m,l) = ', dlm_i(m,l), '; m = ', m, '; l = ', l
715		end do
716
717		do lm = 1, ShapeMap%Shapes(st1)%nRangeFuncs
#	Line 667 | Line 719 | contains
719		m = ShapeMap%Shapes(st1)%RangeFuncMValue(lm)
720		coeff = ShapeMap%Shapes(st1)%RangeFuncCoefficient(lm)
721		function_type = ShapeMap%Shapes(st1)%RangeFunctionType(lm)
722	<
722	>
723		if ((function_type .eq. SH_COS).or.(m.eq.0)) then
724		Phunc = coeff * tm_i(m)
725		dPhuncdX = coeff * dtm_i(m) * dcpidx
726		dPhuncdY = coeff * dtm_i(m) * dcpidy
727		dPhuncdZ = coeff * dtm_i(m) * dcpidz
728	<	dPhuncdUz = coeff * dtm_i(m) * dcpidux
728	>	dPhuncdUx = coeff * dtm_i(m) * dcpidux
729		dPhuncdUy = coeff * dtm_i(m) * dcpiduy
730		dPhuncdUz = coeff * dtm_i(m) * dcpiduz
731		else
#	Line 686 | Line 738 | contains
738		dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
739		endif
740
741	<	s_i = s_i + plm_i(l,m)*Phunc
690	<
691	<	dsidx = dsidx + plm_i(l,m)*dPhuncdX + &
692	<	Phunc * dlm_i(l,m) * dctidx
693	<	dsidy = dsidy + plm_i(l,m)*dPhuncdY + &
694	<	Phunc * dlm_i(l,m) * dctidy
695	<	dsidz = dsidz + plm_i(l,m)*dPhuncdZ + &
696	<	Phunc * dlm_i(l,m) * dctidz
697	<
698	<	dsidux = dsidux + plm_i(l,m)* dPhuncdUx + &
699	<	Phunc * dlm_i(l,m) * dctidux
700	<	dsiduy = dsiduy + plm_i(l,m)* dPhuncdUy + &
701	<	Phunc * dlm_i(l,m) * dctiduy
702	<	dsiduz = dsiduz + plm_i(l,m)* dPhuncdUz + &
703	<	Phunc * dlm_i(l,m) * dctiduz
741	>	s_i = s_i + plm_i(m,l)*Phunc
742
743	+	dsidx = dsidx + plm_i(m,l)*dPhuncdX + &
744	+	Phunc * dlm_i(m,l) * dctidx
745	+	dsidy = dsidy + plm_i(m,l)*dPhuncdY + &
746	+	Phunc * dlm_i(m,l) * dctidy
747	+	dsidz = dsidz + plm_i(m,l)*dPhuncdZ + &
748	+	Phunc * dlm_i(m,l) * dctidz
749	+
750	+	dsidux = dsidux + plm_i(m,l)* dPhuncdUx + &
751	+	Phunc * dlm_i(m,l) * dctidux
752	+	dsiduy = dsiduy + plm_i(m,l)* dPhuncdUy + &
753	+	Phunc * dlm_i(m,l) * dctiduy
754	+	dsiduz = dsiduz + plm_i(m,l)* dPhuncdUz + &
755	+	Phunc * dlm_i(m,l) * dctiduz
756	+
757		end do
758	<
758	>
759		do lm = 1, ShapeMap%Shapes(st1)%nStrengthFuncs
760		l = ShapeMap%Shapes(st1)%StrengthFuncLValue(lm)
761		m = ShapeMap%Shapes(st1)%StrengthFuncMValue(lm)
762		coeff = ShapeMap%Shapes(st1)%StrengthFuncCoefficient(lm)
763		function_type = ShapeMap%Shapes(st1)%StrengthFunctionType(lm)
764	<
764	>
765		if ((function_type .eq. SH_COS).or.(m.eq.0)) then
766		Phunc = coeff * tm_i(m)
767		dPhuncdX = coeff * dtm_i(m) * dcpidx
768		dPhuncdY = coeff * dtm_i(m) * dcpidy
769		dPhuncdZ = coeff * dtm_i(m) * dcpidz
770	<	dPhuncdUz = coeff * dtm_i(m) * dcpidux
770	>	dPhuncdUx = coeff * dtm_i(m) * dcpidux
771		dPhuncdUy = coeff * dtm_i(m) * dcpiduy
772		dPhuncdUz = coeff * dtm_i(m) * dcpiduz
773		else
#	Line 728 | Line 780 | contains
780		dPhuncdUz = coeff*(spi * dum_i(m-1)*dcpiduz + dspiduz *um_i(m-1))
781		endif
782
783	<	eps_i = eps_i + plm_i(l,m)*Phunc
732	<
733	<	depsidx = depsidx + plm_i(l,m)*dPhuncdX + &
734	<	Phunc * dlm_i(l,m) * dctidx
735	<	depsidy = depsidy + plm_i(l,m)*dPhuncdY + &
736	<	Phunc * dlm_i(l,m) * dctidy
737	<	depsidz = depsidz + plm_i(l,m)*dPhuncdZ + &
738	<	Phunc * dlm_i(l,m) * dctidz
739	<
740	<	depsidux = depsidux + plm_i(l,m)* dPhuncdUx + &
741	<	Phunc * dlm_i(l,m) * dctidux
742	<	depsiduy = depsiduy + plm_i(l,m)* dPhuncdUy + &
743	<	Phunc * dlm_i(l,m) * dctiduy
744	<	depsiduz = depsiduz + plm_i(l,m)* dPhuncdUz + &
745	<	Phunc * dlm_i(l,m) * dctiduz
783	>	eps_i = eps_i + plm_i(m,l)*Phunc
784
785	+	depsidx = depsidx + plm_i(m,l)*dPhuncdX + &
786	+	Phunc * dlm_i(m,l) * dctidx
787	+	depsidy = depsidy + plm_i(m,l)*dPhuncdY + &
788	+	Phunc * dlm_i(m,l) * dctidy
789	+	depsidz = depsidz + plm_i(m,l)*dPhuncdZ + &
790	+	Phunc * dlm_i(m,l) * dctidz
791	+
792	+	depsidux = depsidux + plm_i(m,l)* dPhuncdUx + &
793	+	Phunc * dlm_i(m,l) * dctidux
794	+	depsiduy = depsiduy + plm_i(m,l)* dPhuncdUy + &
795	+	Phunc * dlm_i(m,l) * dctiduy
796	+	depsiduz = depsiduz + plm_i(m,l)* dPhuncdUz + &
797	+	Phunc * dlm_i(m,l) * dctiduz
798	+
799		end do
800
801		endif
750	–
751	–	! now do j:
802
803	+	! now do j:
804	+
805		if (ShapeMap%Shapes(st2)%isLJ) then
806		sigma_j = ShapeMap%Shapes(st2)%sigma
807		s_j = ShapeMap%Shapes(st2)%sigma
#	Line 773 | Line 825 | contains
825		depsjduy = 0.0d0
826		depsjduz = 0.0d0
827		else
828	<
828	>
829		#ifdef IS_MPI
830		! rotate the inter-particle separation into the two different
831		! body-fixed coordinate systems:
832		! negative sign because this is the vector from j to i:
833	<
833	>
834		xj = -(A_Col(1,atom2)*d(1) + A_Col(2,atom2)*d(2) + A_Col(3,atom2)*d(3))
835		yj = -(A_Col(4,atom2)*d(1) + A_Col(5,atom2)*d(2) + A_Col(6,atom2)*d(3))
836		zj = -(A_Col(7,atom2)*d(1) + A_Col(8,atom2)*d(2) + A_Col(9,atom2)*d(3))
#	Line 786 | Line 838 | contains
838		! rotate the inter-particle separation into the two different
839		! body-fixed coordinate systems:
840		! negative sign because this is the vector from j to i:
841	<
841	>
842		xj = -(a(1,atom2)*d(1) + a(2,atom2)*d(2) + a(3,atom2)*d(3))
843		yj = -(a(4,atom2)*d(1) + a(5,atom2)*d(2) + a(6,atom2)*d(3))
844		zj = -(a(7,atom2)*d(1) + a(8,atom2)*d(2) + a(9,atom2)*d(3))
845		#endif
846	<
846	>
847	>	xjhat = xj / rij
848	>	yjhat = yj / rij
849	>	zjhat = zj / rij
850		xj2 = xj*xj
851		yj2 = yj*yj
852		zj2 = zj*zj
798	–
799	–	projj = sqrt(xj2 + yj2)
800	–	projj3 = projj*projj*projj
801	–
853		ctj = zj / rij
854	+
855	+	if (ctj .gt. 1.0_dp) ctj = 1.0_dp
856	+	if (ctj .lt. -1.0_dp) ctj = -1.0_dp
857	+
858		dctjdx = - zj * xj / r3
859		dctjdy = - zj * yj / r3
860		dctjdz = 1.0d0 / rij - zj2 / r3
861	<	dctjdux =  yj / rij
862	<	dctjduy = -xj / rij
863	<	dctjduz = 0.0d0
864	<
865	<	cpj = xj / projj
866	<	dcpjdx = 1.0d0 / projj - xj2 / projj3
867	<	dcpjdy = - xj * yj / projj3
868	<	dcpjdz = 0.0d0
869	<	dcpjdux = xj * yj * zj / projj3
870	<	dcpjduy = -zj * (1.0d0 / projj - xj2 / projj3)
871	<	dcpjduz = -yj * (1.0d0 / projj - xj2 / projj3)  - (xj2 * yj / projj3)
872	<
873	<	spj = yj / projj
874	<	dspjdx = - xj * yj / projj3
875	<	dspjdy = 1.0d0 / projj - yj2 / projj3
876	<	dspjdz = 0.0d0
877	<	dspjdux = -zj * (1.0d0 / projj - yj2 / projj3)
878	<	dspjduy = xj * yj * zj / projj3
879	<	dspjduz = xj * (1.0d0 / projj - yi2 / projj3) + (xj * yj2 / projj3)
880	<
881	<	call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigL, &
882	<	ShapeMap%Shapes(st2)%bigM, lmax, plm_j, dlm_j)
883	<
884	<	call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
885	<	CHEBYSHEV_TN, tm_j, dtm_j)
886	<	call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, &
861	>	dctjdux = yj / rij !- (zi * xj2) / r3
862	>	dctjduy = -xj / rij !- (zj * yj2) / r3
863	>	dctjduz = 0.0d0 !zj / rij - (zj2 * zj) / r3
864	>
865	>	! this is an attempt to try to truncate the singularity when
866	>	! sin(theta) is near 0.0:
867	>
868	>	stj2 = 1.0_dp - ctj*ctj
869	>	if (dabs(stj2) .lt. 1.0d-12) then
870	>	projj = sqrt(rij * 1.0d-12)
871	>	dcpjdx = 1.0d0 / projj
872	>	dcpjdy = 0.0d0
873	>	dcpjdux = xj / projj
874	>	dcpjduy = 0.0d0
875	>	dspjdx = 0.0d0
876	>	dspjdy = 1.0d0 / projj
877	>	dspjdux = 0.0d0
878	>	dspjduy = yj / projj
879	>	else
880	>	projj = sqrt(xj2 + yj2)
881	>	projj3 = projj*projj*projj
882	>	dcpjdx = 1.0d0 / projj - xj2 / projj3
883	>	dcpjdy = - xj * yj / projj3
884	>	dcpjdux = xj / projj - (xj2 * xj) / projj3
885	>	dcpjduy = - (xj * yj2) / projj3
886	>	dspjdx = - xj * yj / projj3
887	>	dspjdy = 1.0d0 / projj - yj2 / projj3
888	>	dspjdux = - (yj * xj2) / projj3
889	>	dspjduy = yj / projj - (yj2 * yj) / projj3
890	>	endif
891	>
892	>	cpj = xj / projj
893	>	dcpjdz = 0.0d0
894	>	dcpjduz = 0.0d0
895	>
896	>	spj = yj / projj
897	>	dspjdz = 0.0d0
898	>	dspjduz = 0.0d0
899	>
900	>
901	>	!       write(*,*) 'dcpdu = ' ,dcpidux, dcpiduy, dcpiduz
902	>	!       write(*,*) 'dcpdu = ' ,dcpjdux, dcpjduy, dcpjduz
903	>	call Associated_Legendre(ctj, ShapeMap%Shapes(st2)%bigM, &
904	>	ShapeMap%Shapes(st2)%bigL, LMAX, &
905	>	plm_j, dlm_j)
906	>
907	>	call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
908	>	CHEBYSHEV_TN, tm_j, dtm_j)
909	>	call Orthogonal_Polynomial(cpj, ShapeMap%Shapes(st2)%bigM, MMAX, &
910		CHEBYSHEV_UN, um_j, dum_j)
911	<
911	>
912		sigma_j = 0.0d0
913		s_j = 0.0d0
914		eps_j = 0.0d0
#	Line 864 | Line 942 | contains
942		dPhuncdX = coeff * dtm_j(m) * dcpjdx
943		dPhuncdY = coeff * dtm_j(m) * dcpjdy
944		dPhuncdZ = coeff * dtm_j(m) * dcpjdz
945	<	dPhuncdUz = coeff * dtm_j(m) * dcpjdux
945	>	dPhuncdUx = coeff * dtm_j(m) * dcpjdux
946		dPhuncdUy = coeff * dtm_j(m) * dcpjduy
947		dPhuncdUz = coeff * dtm_j(m) * dcpjduz
948		else
#	Line 876 | Line 954 | contains
954		dPhuncdUy = coeff*(spj * dum_j(m-1)*dcpjduy + dspjduy *um_j(m-1))
955		dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
956		endif
879	–
880	–	sigma_j = sigma_j + plm_j(l,m)*Phunc
881	–
882	–	dsigmajdx = dsigmajdx + plm_j(l,m)*dPhuncdX + &
883	–	Phunc * dlm_j(l,m) * dctjdx
884	–	dsigmajdy = dsigmajdy + plm_j(l,m)*dPhuncdY + &
885	–	Phunc * dlm_j(l,m) * dctjdy
886	–	dsigmajdz = dsigmajdz + plm_j(l,m)*dPhuncdZ + &
887	–	Phunc * dlm_j(l,m) * dctjdz
888	–
889	–	dsigmajdux = dsigmajdux + plm_j(l,m)* dPhuncdUx + &
890	–	Phunc * dlm_j(l,m) * dctjdux
891	–	dsigmajduy = dsigmajduy + plm_j(l,m)* dPhuncdUy + &
892	–	Phunc * dlm_j(l,m) * dctjduy
893	–	dsigmajduz = dsigmajduz + plm_j(l,m)* dPhuncdUz + &
894	–	Phunc * dlm_j(l,m) * dctjduz
957
958	+	sigma_j = sigma_j + plm_j(m,l)*Phunc
959	+
960	+	dsigmajdx = dsigmajdx + plm_j(m,l)*dPhuncdX + &
961	+	Phunc * dlm_j(m,l) * dctjdx
962	+	dsigmajdy = dsigmajdy + plm_j(m,l)*dPhuncdY + &
963	+	Phunc * dlm_j(m,l) * dctjdy
964	+	dsigmajdz = dsigmajdz + plm_j(m,l)*dPhuncdZ + &
965	+	Phunc * dlm_j(m,l) * dctjdz
966	+
967	+	dsigmajdux = dsigmajdux + plm_j(m,l)* dPhuncdUx + &
968	+	Phunc * dlm_j(m,l) * dctjdux
969	+	dsigmajduy = dsigmajduy + plm_j(m,l)* dPhuncdUy + &
970	+	Phunc * dlm_j(m,l) * dctjduy
971	+	dsigmajduz = dsigmajduz + plm_j(m,l)* dPhuncdUz + &
972	+	Phunc * dlm_j(m,l) * dctjduz
973	+
974		end do
975
976		do lm = 1, ShapeMap%Shapes(st2)%nRangeFuncs
#	Line 906 | Line 984 | contains
984		dPhuncdX = coeff * dtm_j(m) * dcpjdx
985		dPhuncdY = coeff * dtm_j(m) * dcpjdy
986		dPhuncdZ = coeff * dtm_j(m) * dcpjdz
987	<	dPhuncdUz = coeff * dtm_j(m) * dcpjdux
987	>	dPhuncdUx = coeff * dtm_j(m) * dcpjdux
988		dPhuncdUy = coeff * dtm_j(m) * dcpjduy
989		dPhuncdUz = coeff * dtm_j(m) * dcpjduz
990		else
#	Line 919 | Line 997 | contains
997		dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
998		endif
999
1000	<	s_j = s_j + plm_j(l,m)*Phunc
923	<
924	<	dsjdx = dsjdx + plm_j(l,m)*dPhuncdX + &
925	<	Phunc * dlm_j(l,m) * dctjdx
926	<	dsjdy = dsjdy + plm_j(l,m)*dPhuncdY + &
927	<	Phunc * dlm_j(l,m) * dctjdy
928	<	dsjdz = dsjdz + plm_j(l,m)*dPhuncdZ + &
929	<	Phunc * dlm_j(l,m) * dctjdz
930	<
931	<	dsjdux = dsjdux + plm_j(l,m)* dPhuncdUx + &
932	<	Phunc * dlm_j(l,m) * dctjdux
933	<	dsjduy = dsjduy + plm_j(l,m)* dPhuncdUy + &
934	<	Phunc * dlm_j(l,m) * dctjduy
935	<	dsjduz = dsjduz + plm_j(l,m)* dPhuncdUz + &
936	<	Phunc * dlm_j(l,m) * dctjduz
1000	>	s_j = s_j + plm_j(m,l)*Phunc
1001
1002	+	dsjdx = dsjdx + plm_j(m,l)*dPhuncdX + &
1003	+	Phunc * dlm_j(m,l) * dctjdx
1004	+	dsjdy = dsjdy + plm_j(m,l)*dPhuncdY + &
1005	+	Phunc * dlm_j(m,l) * dctjdy
1006	+	dsjdz = dsjdz + plm_j(m,l)*dPhuncdZ + &
1007	+	Phunc * dlm_j(m,l) * dctjdz
1008	+
1009	+	dsjdux = dsjdux + plm_j(m,l)* dPhuncdUx + &
1010	+	Phunc * dlm_j(m,l) * dctjdux
1011	+	dsjduy = dsjduy + plm_j(m,l)* dPhuncdUy + &
1012	+	Phunc * dlm_j(m,l) * dctjduy
1013	+	dsjduz = dsjduz + plm_j(m,l)* dPhuncdUz + &
1014	+	Phunc * dlm_j(m,l) * dctjduz
1015	+
1016		end do
1017
1018		do lm = 1, ShapeMap%Shapes(st2)%nStrengthFuncs
#	Line 961 | Line 1039 | contains
1039		dPhuncdUz = coeff*(spj * dum_j(m-1)*dcpjduz + dspjduz *um_j(m-1))
1040		endif
1041
1042	<	eps_j = eps_j + plm_j(l,m)*Phunc
965	<
966	<	depsjdx = depsjdx + plm_j(l,m)*dPhuncdX + &
967	<	Phunc * dlm_j(l,m) * dctjdx
968	<	depsjdy = depsjdy + plm_j(l,m)*dPhuncdY + &
969	<	Phunc * dlm_j(l,m) * dctjdy
970	<	depsjdz = depsjdz + plm_j(l,m)*dPhuncdZ + &
971	<	Phunc * dlm_j(l,m) * dctjdz
972	<
973	<	depsjdux = depsjdux + plm_j(l,m)* dPhuncdUx + &
974	<	Phunc * dlm_j(l,m) * dctjdux
975	<	depsjduy = depsjduy + plm_j(l,m)* dPhuncdUy + &
976	<	Phunc * dlm_j(l,m) * dctjduy
977	<	depsjduz = depsjduz + plm_j(l,m)* dPhuncdUz + &
978	<	Phunc * dlm_j(l,m) * dctjduz
1042	>	!          write(*,*) 'l,m = ', l, m, coeff, dPhuncdUx, dPhuncdUy, dPhuncdUz
1043
1044	+	eps_j = eps_j + plm_j(m,l)*Phunc
1045	+
1046	+	depsjdx = depsjdx + plm_j(m,l)*dPhuncdX + &
1047	+	Phunc * dlm_j(m,l) * dctjdx
1048	+	depsjdy = depsjdy + plm_j(m,l)*dPhuncdY + &
1049	+	Phunc * dlm_j(m,l) * dctjdy
1050	+	depsjdz = depsjdz + plm_j(m,l)*dPhuncdZ + &
1051	+	Phunc * dlm_j(m,l) * dctjdz
1052	+
1053	+	depsjdux = depsjdux + plm_j(m,l)* dPhuncdUx + &
1054	+	Phunc * dlm_j(m,l) * dctjdux
1055	+	depsjduy = depsjduy + plm_j(m,l)* dPhuncdUy + &
1056	+	Phunc * dlm_j(m,l) * dctjduy
1057	+	depsjduz = depsjduz + plm_j(m,l)* dPhuncdUz + &
1058	+	Phunc * dlm_j(m,l) * dctjduz
1059	+
1060		end do
1061
1062		endif
#	Line 984 | Line 1064 | contains
1064		! phew, now let's assemble the potential energy:
1065
1066		sigma = 0.5*(sigma_i + sigma_j)
1067	<
1067	>	!    write(*,*) sigma_i, ' = sigma_i; ', sigma_j, ' = sigma_j'
1068		dsigmadxi = 0.5*dsigmaidx
1069		dsigmadyi = 0.5*dsigmaidy
1070		dsigmadzi = 0.5*dsigmaidz
#	Line 1016 | Line 1096 | contains
1096		dsduzj = 0.5*dsjduz
1097
1098		eps = sqrt(eps_i * eps_j)
1099	<
1099	>	!!$    write(*,*) 'dsidu = ', dsidux, dsiduy, dsiduz
1100	>	!!$    write(*,*) 'dsigidu = ', dsigmaidux, dsigmaiduy, dsigmaiduz
1101	>	!!$    write(*,*) sigma_j, ' is sigma j; ', s_j, ' is s j; ', eps_j, ' is eps j'
1102		depsdxi = eps_j * depsidx / (2.0d0 * eps)
1103		depsdyi = eps_j * depsidy / (2.0d0 * eps)
1104		depsdzi = eps_j * depsidz / (2.0d0 * eps)
#	Line 1030 | Line 1112 | contains
1112		depsduxj = eps_i * depsjdux / (2.0d0 * eps)
1113		depsduyj = eps_i * depsjduy / (2.0d0 * eps)
1114		depsduzj = eps_i * depsjduz / (2.0d0 * eps)
1115	<
1115	>
1116	>	!!$    write(*,*) 'depsidu = ', depsidux, depsiduy, depsiduz
1117	>
1118	>	!!$    write(*,*) 'depsjdu = ', depsjdux, depsjduy, depsjduz
1119	>	!!$    write(*,*) 'depsduj = ', depsduxj, depsduyj, depsduzj
1120	>	!!$
1121	>	!!$    write(*,*) 's, sig, eps = ', s, sigma, eps
1122	>
1123		rtdenom = rij-sigma+s
1124		rt = s / rtdenom
1125
1126	<	drtdxi = (dsdxi + rt * (drdxi - dsigmadxi + dsdxi)) / rtdenom
1127	<	drtdyi = (dsdyi + rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1128	<	drtdzi = (dsdzi + rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
1129	<	drtduxi = (dsduxi + rt * (drduxi - dsigmaduxi + dsduxi)) / rtdenom
1130	<	drtduyi = (dsduyi + rt * (drduyi - dsigmaduyi + dsduyi)) / rtdenom
1131	<	drtduzi = (dsduzi + rt * (drduzi - dsigmaduzi + dsduzi)) / rtdenom
1132	<	drtdxj = (dsdxj + rt * (drdxj - dsigmadxj + dsdxj)) / rtdenom
1133	<	drtdyj = (dsdyj + rt * (drdyj - dsigmadyj + dsdyj)) / rtdenom
1134	<	drtdzj = (dsdzj + rt * (drdzj - dsigmadzj + dsdzj)) / rtdenom
1135	<	drtduxj = (dsduxj + rt * (drduxj - dsigmaduxj + dsduxj)) / rtdenom
1136	<	drtduyj = (dsduyj + rt * (drduyj - dsigmaduyj + dsduyj)) / rtdenom
1137	<	drtduzj = (dsduzj + rt * (drduzj - dsigmaduzj + dsduzj)) / rtdenom
1138	<
1126	>	drtdxi = (dsdxi - rt * (drdxi - dsigmadxi + dsdxi)) / rtdenom
1127	>	drtdyi = (dsdyi - rt * (drdyi - dsigmadyi + dsdyi)) / rtdenom
1128	>	drtdzi = (dsdzi - rt * (drdzi - dsigmadzi + dsdzi)) / rtdenom
1129	>	drtduxi = (dsduxi - rt * (drduxi - dsigmaduxi + dsduxi)) / rtdenom
1130	>	drtduyi = (dsduyi - rt * (drduyi - dsigmaduyi + dsduyi)) / rtdenom
1131	>	drtduzi = (dsduzi - rt * (drduzi - dsigmaduzi + dsduzi)) / rtdenom
1132	>	drtdxj = (dsdxj - rt * (drdxj - dsigmadxj + dsdxj)) / rtdenom
1133	>	drtdyj = (dsdyj - rt * (drdyj - dsigmadyj + dsdyj)) / rtdenom
1134	>	drtdzj = (dsdzj - rt * (drdzj - dsigmadzj + dsdzj)) / rtdenom
1135	>	drtduxj = (dsduxj - rt * (drduxj - dsigmaduxj + dsduxj)) / rtdenom
1136	>	drtduyj = (dsduyj - rt * (drduyj - dsigmaduyj + dsduyj)) / rtdenom
1137	>	drtduzj = (dsduzj - rt * (drduzj - dsigmaduzj + dsduzj)) / rtdenom
1138	>
1139	>	!!$    write(*,*) 'drtd_i = ', drtdxi, drtdyi, drtdzi
1140	>	!!$    write(*,*) 'drtdu_j = ', drtduxj, drtduyj, drtduzj
1141	>
1142		rt2 = rt*rt
1143		rt3 = rt2*rt
1144		rt5 = rt2*rt3
#	Line 1055 | Line 1147 | contains
1147		rt12 = rt6*rt6
1148		rt126 = rt12 - rt6
1149
1150	+	pot_temp = 4.0d0 * eps * rt126
1151	+
1152	+	vpair = vpair + pot_temp
1153		if (do_pot) then
1154		#ifdef IS_MPI
1155	<	pot_row(atom1) = pot_row(atom1) + 2.0d0*eps*rt126*sw
1156	<	pot_col(atom2) = pot_col(atom2) + 2.0d0*eps*rt126*sw
1155	>	pot_row(SHAPE_POT,atom1) = pot_row(SHAPE_POT,atom1) + 0.5d0*pot_temp*sw
1156	>	pot_col(SHAPE_POT,atom2) = pot_col(SHAPE_POT,atom2) + 0.5d0*pot_temp*sw
1157		#else
1158	<	pot = pot + 4.0d0*eps*rt126*sw
1158	>	pot = pot + pot_temp*sw
1159		#endif
1160		endif
1161	<
1161	>
1162	>	!!$    write(*,*) 'drtdu, depsdu = ', drtduxi, depsduxi
1163	>
1164		dvdxi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdxi + 4.0d0*depsdxi*rt126
1165		dvdyi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdyi + 4.0d0*depsdyi*rt126
1166		dvdzi = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtdzi + 4.0d0*depsdzi*rt126
#	Line 1077 | Line 1174 | contains
1174		dvduxj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduxj + 4.0d0*depsduxj*rt126
1175		dvduyj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduyj + 4.0d0*depsduyj*rt126
1176		dvduzj = 24.0d0*eps*(2.0d0*rt11 - rt5)*drtduzj + 4.0d0*depsduzj*rt126
1177	<
1177	>	!!$    write(*,*) 'drtduxi = ', drtduxi, ' depsduxi = ', depsduxi
1178		! do the torques first since they are easy:
1179		! remember that these are still in the body fixed axes
1180
1181	<	txi = dvduxi * sw
1182	<	tyi = dvduyi * sw
1183	<	tzi = dvduzi * sw
1181	>	txi = 0.0d0
1182	>	tyi = 0.0d0
1183	>	tzi = 0.0d0
1184
1185	<	txj = dvduxj * sw
1186	<	tyj = dvduyj * sw
1187	<	tzj = dvduzj * sw
1185	>	txj = 0.0d0
1186	>	tyj = 0.0d0
1187	>	tzj = 0.0d0
1188
1189	+	txi = (dvduyi - dvduzi) * sw
1190	+	tyi = (dvduzi - dvduxi) * sw
1191	+	tzi = (dvduxi - dvduyi) * sw
1192	+
1193	+	txj = (dvduyj - dvduzj) * sw
1194	+	tyj = (dvduzj - dvduxj) * sw
1195	+	tzj = (dvduxj - dvduyj) * sw
1196	+
1197	+	!!$    txi = dvduxi * sw
1198	+	!!$    tyi = dvduyi * sw
1199	+	!!$    tzi = dvduzi * sw
1200	+	!!$
1201	+	!!$    txj = dvduxj * sw
1202	+	!!$    tyj = dvduyj * sw
1203	+	!!$    tzj = dvduzj * sw
1204	+
1205	+	write(*,*) 't1 = ', txi, tyi, tzi
1206	+	write(*,*) 't2 = ', txj, tyj, tzj
1207	+
1208		! go back to lab frame using transpose of rotation matrix:
1209	<
1209	>
1210		#ifdef IS_MPI
1211		t_Row(1,atom1) = t_Row(1,atom1) + a_Row(1,atom1)*txi + &
1212		a_Row(4,atom1)*tyi + a_Row(7,atom1)*tzi
#	Line 1098 | Line 1214 | contains
1214		a_Row(5,atom1)*tyi + a_Row(8,atom1)*tzi
1215		t_Row(3,atom1) = t_Row(3,atom1) + a_Row(3,atom1)*txi + &
1216		a_Row(6,atom1)*tyi + a_Row(9,atom1)*tzi
1217	<
1217	>
1218		t_Col(1,atom2) = t_Col(1,atom2) + a_Col(1,atom2)*txj + &
1219		a_Col(4,atom2)*tyj + a_Col(7,atom2)*tzj
1220		t_Col(2,atom2) = t_Col(2,atom2) + a_Col(2,atom2)*txj + &
1221	<	a_Col(5,atom2)*tyj + a_Col(8,atom2)*tzj
1221	>	a_Col(5,atom2)*tyj + a_Col(8,atom2)*tzj
1222		t_Col(3,atom2) = t_Col(3,atom2) + a_Col(3,atom2)*txj + &
1223		a_Col(6,atom2)*tyj + a_Col(9,atom2)*tzj
1224		#else
1225		t(1,atom1) = t(1,atom1) + a(1,atom1)*txi + a(4,atom1)*tyi + a(7,atom1)*tzi
1226		t(2,atom1) = t(2,atom1) + a(2,atom1)*txi + a(5,atom1)*tyi + a(8,atom1)*tzi
1227		t(3,atom1) = t(3,atom1) + a(3,atom1)*txi + a(6,atom1)*tyi + a(9,atom1)*tzi
1228	<
1228	>
1229		t(1,atom2) = t(1,atom2) + a(1,atom2)*txj + a(4,atom2)*tyj + a(7,atom2)*tzj
1230		t(2,atom2) = t(2,atom2) + a(2,atom2)*txj + a(5,atom2)*tyj + a(8,atom2)*tzj
1231		t(3,atom2) = t(3,atom2) + a(3,atom2)*txj + a(6,atom2)*tyj + a(9,atom2)*tzj
1232	+
1233		#endif
1234		! Now, on to the forces:
1235	<
1235	>
1236		! first rotate the i terms back into the lab frame:
1120	–
1121	–	fxi = dvdxi * sw
1122	–	fyi = dvdyi * sw
1123	–	fzi = dvdzi * sw
1237
1238	<	fxj = dvdxj * sw
1239	<	fyj = dvdyj * sw
1240	<	fzj = dvdzj * sw
1238	>	fxi = -dvdxi * sw
1239	>	fyi = -dvdyi * sw
1240	>	fzi = -dvdzi * sw
1241
1242	+	fxj = -dvdxj * sw
1243	+	fyj = -dvdyj * sw
1244	+	fzj = -dvdzj * sw
1245	+
1246	+
1247		#ifdef IS_MPI
1248		fxii = a_Row(1,atom1)*fxi + a_Row(4,atom1)*fyi + a_Row(7,atom1)*fzi
1249		fyii = a_Row(2,atom1)*fxi + a_Row(5,atom1)*fyi + a_Row(8,atom1)*fzi
#	Line 1138 | Line 1256 | contains
1256		fxii = a(1,atom1)*fxi + a(4,atom1)*fyi + a(7,atom1)*fzi
1257		fyii = a(2,atom1)*fxi + a(5,atom1)*fyi + a(8,atom1)*fzi
1258		fzii = a(3,atom1)*fxi + a(6,atom1)*fyi + a(9,atom1)*fzi
1259	<
1259	>
1260		fxjj = a(1,atom2)*fxj + a(4,atom2)*fyj + a(7,atom2)*fzj
1261		fyjj = a(2,atom2)*fxj + a(5,atom2)*fyj + a(8,atom2)*fzj
1262		fzjj = a(3,atom2)*fxj + a(6,atom2)*fyj + a(9,atom2)*fzj
#	Line 1147 | Line 1265 | contains
1265		fxij = -fxii
1266		fyij = -fyii
1267		fzij = -fzii
1268	<
1268	>
1269		fxji = -fxjj
1270		fyji = -fyjj
1271		fzji = -fzjj
1272
1273	<	fxradial = fxii + fxji
1274	<	fyradial = fyii + fyji
1275	<	fzradial = fzii + fzji
1276	<
1273	>	fxradial = (fxii + fxji)
1274	>	fyradial = (fyii + fyji)
1275	>	fzradial = (fzii + fzji)
1276	>	!!$    write(*,*) fxradial, ' is fxrad; ', fyradial, ' is fyrad; ', fzradial, 'is fzrad'
1277		#ifdef IS_MPI
1278		f_Row(1,atom1) = f_Row(1,atom1) + fxradial
1279		f_Row(2,atom1) = f_Row(2,atom1) + fyradial
1280		f_Row(3,atom1) = f_Row(3,atom1) + fzradial
1281	<
1281	>
1282		f_Col(1,atom2) = f_Col(1,atom2) - fxradial
1283		f_Col(2,atom2) = f_Col(2,atom2) - fyradial
1284		f_Col(3,atom2) = f_Col(3,atom2) - fzradial
#	Line 1168 | Line 1286 | contains
1286		f(1,atom1) = f(1,atom1) + fxradial
1287		f(2,atom1) = f(2,atom1) + fyradial
1288		f(3,atom1) = f(3,atom1) + fzradial
1289	<
1289	>
1290		f(1,atom2) = f(1,atom2) - fxradial
1291		f(2,atom2) = f(2,atom2) - fyradial
1292		f(3,atom2) = f(3,atom2) - fzradial
#	Line 1181 | Line 1299 | contains
1299		id1 = atom1
1300		id2 = atom2
1301		#endif
1302	<
1302	>
1303		if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1304	<
1304	>
1305		fpair(1) = fpair(1) + fxradial
1306		fpair(2) = fpair(2) + fyradial
1307		fpair(3) = fpair(3) + fzradial
1308	<
1308	>
1309		endif
1310	<
1310	>
1311		end subroutine do_shape_pair
1312	<
1313	<	SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)
1314	<
1312	>
1313	>	SUBROUTINE Associated_Legendre(x, l, m, lmax, plm, dlm)
1314	>
1315		! Purpose: Compute the associated Legendre functions
1316		!          Plm(x) and their derivatives Plm'(x)
1317		! Input :  x  --- Argument of Plm(x)
#	Line 1209 | Line 1327 | contains
1327		!
1328		! The original Fortran77 codes can be found here:
1329		! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1330	<
1331	<	real (kind=8), intent(in) :: x
1330	>
1331	>	real (kind=dp), intent(in) :: x
1332		integer, intent(in) :: l, m, lmax
1333	<	real (kind=8), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
1333	>	real (kind=dp), dimension(0:lmax,0:m), intent(out) :: PLM, DLM
1334		integer :: i, j, ls
1335	<	real (kind=8) :: xq, xs
1335	>	real (kind=dp) :: xq, xs
1336
1337		! zero out both arrays:
1338		DO I = 0, m
1339		DO J = 0, l
1340	<	PLM(J,I) = 0.0D0
1341	<	DLM(J,I) = 0.0D0
1340	>	PLM(J,I) = 0.0_dp
1341	>	DLM(J,I) = 0.0_dp
1342		end DO
1343		end DO
1344
1345		! start with 0,0:
1346		PLM(0,0) = 1.0D0
1347	<
1347	>
1348		! x = +/- 1 functions are easy:
1349		IF (abs(X).EQ.1.0D0) THEN
1350		DO I = 1, m
#	Line 1253 | Line 1371 | contains
1371		DO I = 1, l
1372		PLM(I, I) = -LS*(2.0D0*I-1.0D0)*XQ*PLM(I-1, I-1)
1373		enddo
1374	<
1374	>
1375		DO I = 0, l
1376		PLM(I, I+1)=(2.0D0*I+1.0D0)*X*PLM(I, I)
1377		enddo
1378	<
1378	>
1379		DO I = 0, l
1380		DO J = I+2, m
1381		PLM(I, J)=((2.0D0*J-1.0D0)*X*PLM(I,J-1) - &
1382		(I+J-1.0D0)*PLM(I,J-2))/(J-I)
1383		end DO
1384		end DO
1385	<
1385	>
1386		DLM(0, 0)=0.0D0
1269	–
1387		DO J = 1, m
1388		DLM(0, J)=LS*J*(PLM(0,J-1)-X*PLM(0,J))/XS
1389		end DO
1390	<
1390	>
1391		DO I = 1, l
1392		DO J = I, m
1393		DLM(I,J) = LS*I*X*PLM(I, J)/XS + (J+I)*(J-I+1.0D0)/XQ*PLM(I-1, J)
1394		end DO
1395		end DO
1396	<
1396	>
1397		RETURN
1398		END SUBROUTINE Associated_Legendre
1399
1400
1401	<	subroutine Orthogonal_Polynomial(x, m, function_type, pl, dpl)
1402	<
1401	>	subroutine Orthogonal_Polynomial(x, m, mmax, function_type, pl, dpl)
1402	>
1403		! Purpose: Compute orthogonal polynomials: Tn(x) or Un(x),
1404		!          or Ln(x) or Hn(x), and their derivatives
1405		! Input :  function_type --- Function code
#	Line 1301 | Line 1418 | contains
1418		!
1419		! The original Fortran77 codes can be found here:
1420		! http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
1421	<
1421	>
1422		real(kind=8), intent(in) :: x
1423	<	integer, intent(in):: m
1423	>	integer, intent(in):: m, mmax
1424		integer, intent(in):: function_type
1425	<	real(kind=8), dimension(0:m), intent(inout) :: pl, dpl
1426	<
1425	>	real(kind=8), dimension(0:mmax), intent(inout) :: pl, dpl
1426	>
1427		real(kind=8) :: a, b, c, y0, y1, dy0, dy1, yn, dyn
1428		integer :: k
1429
#	Line 1349 | Line 1466 | contains
1466		DY0 = DY1
1467		DY1 = DYN
1468		end DO
1469	+
1470	+
1471		RETURN
1472	<
1472	>
1473		end subroutine Orthogonal_Polynomial
1355	–
1356	–	end module shapes
1474
1475	<	subroutine makeShape(nContactFuncs, ContactFuncLValue, &
1476	<	ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
1360	<	nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1361	<	RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1362	<	StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1363	<	myAtid, status)
1475	>	subroutine deallocateShapes(this)
1476	>	type(Shape), pointer :: this
1477
1478	<	use definitions
1479	<	use shapes, only: newShapeType
1480	<
1481	<	integer :: nContactFuncs
1369	<	integer :: nRangeFuncs
1370	<	integer :: nStrengthFuncs
1371	<	integer :: status
1372	<	integer :: myAtid
1373	<
1374	<	integer, dimension(nContactFuncs) :: ContactFuncLValue
1375	<	integer, dimension(nContactFuncs) :: ContactFuncMValue
1376	<	integer, dimension(nContactFuncs) :: ContactFunctionType
1377	<	real(kind=dp), dimension(nContactFuncs) :: ContactFuncCoefficient
1378	<	integer, dimension(nRangeFuncs) :: RangeFuncLValue
1379	<	integer, dimension(nRangeFuncs) :: RangeFuncMValue
1380	<	integer, dimension(nRangeFuncs) :: RangeFunctionType
1381	<	real(kind=dp), dimension(nRangeFuncs) :: RangeFuncCoefficient
1382	<	integer, dimension(nStrengthFuncs) :: StrengthFuncLValue
1383	<	integer, dimension(nStrengthFuncs) :: StrengthFuncMValue
1384	<	integer, dimension(nStrengthFuncs) :: StrengthFunctionType
1385	<	real(kind=dp), dimension(nStrengthFuncs) :: StrengthFuncCoefficient
1386	<
1387	<	call newShapeType(nContactFuncs, ContactFuncLValue, &
1388	<	ContactFuncMValue, ContactFunctionType, ContactFuncCoefficient, &
1389	<	nRangeFuncs, RangeFuncLValue, RangeFuncMValue, RangeFunctionType, &
1390	<	RangeFuncCoefficient, nStrengthFuncs, StrengthFuncLValue, &
1391	<	StrengthFuncMValue, StrengthFunctionType, StrengthFuncCoefficient, &
1392	<	myAtid, status)
1478	>	if (associated( this%ContactFuncLValue)) then
1479	>	deallocate(this%ContactFuncLValue)
1480	>	this%ContactFuncLValue => null()
1481	>	end if
1482
1483	<	return
1484	<	end subroutine makeShape
1483	>	if (associated( this%ContactFuncMValue)) then
1484	>	deallocate( this%ContactFuncMValue)
1485	>	this%ContactFuncMValue => null()
1486	>	end if
1487	>	if (associated( this%ContactFunctionType)) then
1488	>	deallocate(this%ContactFunctionType)
1489	>	this%ContactFunctionType => null()
1490	>	end if
1491	>
1492	>	if (associated( this%ContactFuncCoefficient)) then
1493	>	deallocate(this%ContactFuncCoefficient)
1494	>	this%ContactFuncCoefficient => null()
1495	>	end if
1496	>
1497	>	if (associated( this%RangeFuncLValue)) then
1498	>	deallocate(this%RangeFuncLValue)
1499	>	this%RangeFuncLValue => null()
1500	>	end if
1501	>	if (associated( this%RangeFuncMValue)) then
1502	>	deallocate( this%RangeFuncMValue)
1503	>	this%RangeFuncMValue => null()
1504	>	end if
1505	>
1506	>	if (associated( this%RangeFunctionType)) then
1507	>	deallocate( this%RangeFunctionType)
1508	>	this%RangeFunctionType => null()
1509	>	end if
1510	>	if (associated( this%RangeFuncCoefficient)) then
1511	>	deallocate(this%RangeFuncCoefficient)
1512	>	this%RangeFuncCoefficient => null()
1513	>	end if
1514	>
1515	>	if (associated( this%StrengthFuncLValue)) then
1516	>	deallocate(this%StrengthFuncLValue)
1517	>	this%StrengthFuncLValue => null()
1518	>	end if
1519	>
1520	>	if (associated( this%StrengthFuncMValue )) then
1521	>	deallocate(this%StrengthFuncMValue)
1522	>	this%StrengthFuncMValue => null()
1523	>	end if
1524	>
1525	>	if(associated( this%StrengthFunctionType)) then
1526	>	deallocate(this%StrengthFunctionType)
1527	>	this%StrengthFunctionType => null()
1528	>	end if
1529	>	if (associated( this%StrengthFuncCoefficient )) then
1530	>	deallocate(this%StrengthFuncCoefficient)
1531	>	this%StrengthFuncCoefficient => null()
1532	>	end if
1533	>	end subroutine deallocateShapes
1534	>
1535	>	subroutine destroyShapeTypes
1536	>	integer :: i
1537	>	type(Shape), pointer :: thisShape
1538	>
1539	>	! First walk through and kill the shape
1540	>	do i = 1,ShapeMap%n_shapes
1541	>	thisShape => ShapeMap%Shapes(i)
1542	>	call deallocateShapes(thisShape)
1543	>	end do
1544	>
1545	>	! set shape map to starting values
1546	>	ShapeMap%n_shapes = 0
1547	>	ShapeMap%currentShape = 0
1548	>
1549	>	if (associated(ShapeMap%Shapes)) then
1550	>	deallocate(ShapeMap%Shapes)
1551	>	ShapeMap%Shapes => null()
1552	>	end if
1553	>
1554	>	if (associated(ShapeMap%atidToShape)) then
1555	>	deallocate(ShapeMap%atidToShape)
1556	>	ShapeMap%atidToShape => null()
1557	>	end if
1558	>
1559	>
1560	>	end subroutine destroyShapeTypes
1561	>
1562	>
1563	>	end module shapes

Diff Legend

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