[ViewVC] Diff of: group/trunk/OOPSE-2.0/src/UseTheForce/doForces.F90

Comparing trunk/OOPSE-2.0/src/UseTheForce/doForces.F90 (file contents):
Revision 2284 by gezelter, Wed Sep 7 19:18:54 2005 UTC vs.
Revision 2356 by chuckv, Wed Oct 12 19:55:26 2005 UTC

#	Line 45 \| Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.37 2005-09-07 19:18:54 gezelter Exp $, $Date: 2005-09-07 19:18:54 $, $Name: not supported by cvs2svn $, $Revision: 1.37 $
48	>	!! @version $Id: doForces.F90,v 1.55 2005-10-12 19:55:26 chuckv Exp $, $Date: 2005-10-12 19:55:26 $, $Name: not supported by cvs2svn $, $Revision: 1.55 $
49
50
51		module doForces
#	Line 58 \| Line 58 \| module doForces
58		use lj
59		use sticky
60		use electrostatic_module
61	<	use reaction_field
61	>	use reaction_field_module
62		use gb_pair
63		use shapes
64		use vector_class
#	Line 75 \| Line 75 \| module doForces
75		#include "UseTheForce/fSwitchingFunction.h"
76		#include "UseTheForce/fCutoffPolicy.h"
77		#include "UseTheForce/DarkSide/fInteractionMap.h"
78	+	#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
79
80
81		INTEGER, PARAMETER:: PREPAIR_LOOP = 1
#	Line 85 \| Line 86 \| module doForces
86		logical, save :: haveSaneForceField = .false.
87		logical, save :: haveInteractionHash = .false.
88		logical, save :: haveGtypeCutoffMap = .false.
89	+	logical, save :: haveDefaultCutoffs = .false.
90		logical, save :: haveRlist = .false.
91
92		logical, save :: FF_uses_DirectionalAtoms
93		logical, save :: FF_uses_Dipoles
94		logical, save :: FF_uses_GayBerne
95		logical, save :: FF_uses_EAM
94	–	logical, save :: FF_uses_RF
96
97		logical, save :: SIM_uses_DirectionalAtoms
98		logical, save :: SIM_uses_EAM
98	–	logical, save :: SIM_uses_RF
99		logical, save :: SIM_requires_postpair_calc
100		logical, save :: SIM_requires_prepair_calc
101		logical, save :: SIM_uses_PBC
102
103	<	integer, save :: corrMethod
103	>	integer, save :: electrostaticSummationMethod
104
105		public :: init_FF
106		public :: setDefaultCutoffs
#	Line 124 \| Line 124 \| module doForces
124		! Bit hash to determine pair-pair interactions.
125		integer, dimension(:,:), allocatable :: InteractionHash
126		real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
127	<	real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
128	<	integer, dimension(:), allocatable :: groupToGtype
129	<	real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
127	>	real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
128	>	real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
129	>
130	>	integer, dimension(:), allocatable, target :: groupToGtypeRow
131	>	integer, dimension(:), pointer :: groupToGtypeCol => null()
132	>
133	>	real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
134	>	real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
135		type ::gtypeCutoffs
136		real(kind=dp) :: rcut
137		real(kind=dp) :: rcutsq
#	Line 136 \| Line 141 \| module doForces
141
142		integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
143		real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
144	+	real(kind=dp),save :: listSkin
145
146		contains
147
#	Line 179 \| Line 185 \| contains
185
186		if (.not. allocated(InteractionHash)) then
187		allocate(InteractionHash(nAtypes,nAtypes))
188	+	else
189	+	deallocate(InteractionHash)
190	+	allocate(InteractionHash(nAtypes,nAtypes))
191		endif
192
193		if (.not. allocated(atypeMaxCutoff)) then
194		allocate(atypeMaxCutoff(nAtypes))
195	+	else
196	+	deallocate(atypeMaxCutoff)
197	+	allocate(atypeMaxCutoff(nAtypes))
198		endif
199
200		do i = 1, nAtypes
#	Line 256 \| Line 268 \| contains
268		logical :: i_is_GB
269		logical :: i_is_EAM
270		logical :: i_is_Shape
271	+	logical :: GtypeFound
272
273		integer :: myStatus, nAtypes, i, j, istart, iend, jstart, jend
274	<	integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
275	<	real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
274	>	integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
275	>	integer :: nGroupsInRow
276	>	integer :: nGroupsInCol
277	>	integer :: nGroupTypesRow,nGroupTypesCol
278	>	real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
279		real(kind=dp) :: biggestAtypeCutoff
280
281		stat = 0
#	Line 271 \| Line 287 \| contains
287		return
288		endif
289		endif
290	<
290	>	#ifdef IS_MPI
291	>	nGroupsInRow = getNgroupsInRow(plan_group_row)
292	>	nGroupsInCol = getNgroupsInCol(plan_group_col)
293	>	#endif
294		nAtypes = getSize(atypes)
295	<
295	>	! Set all of the initial cutoffs to zero.
296	>	atypeMaxCutoff = 0.0_dp
297		do i = 1, nAtypes
298		if (SimHasAtype(i)) then
299		call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
#	Line 284 \| Line 304 \| contains
304		call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
305		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
306
307	<	atypeMaxCutoff(i) = 0.0_dp
308	<	if (i_is_LJ) then
309	<	thisRcut = getSigma(i) * 2.5_dp
310	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
311	<	endif
312	<	if (i_is_Elect) then
313	<	thisRcut = defaultRcut
314	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
307	>
308	>	if (haveDefaultCutoffs) then
309	>	atypeMaxCutoff(i) = defaultRcut
310	>	else
311	>	if (i_is_LJ) then
312	>	thisRcut = getSigma(i) * 2.5_dp
313	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
314	>	endif
315	>	if (i_is_Elect) then
316	>	thisRcut = defaultRcut
317	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
318	>	endif
319	>	if (i_is_Sticky) then
320	>	thisRcut = getStickyCut(i)
321	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
322	>	endif
323	>	if (i_is_StickyP) then
324	>	thisRcut = getStickyPowerCut(i)
325	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
326	>	endif
327	>	if (i_is_GB) then
328	>	thisRcut = getGayBerneCut(i)
329	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
330	>	endif
331	>	if (i_is_EAM) then
332	>	thisRcut = getEAMCut(i)
333	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
334	>	endif
335	>	if (i_is_Shape) then
336	>	thisRcut = getShapeCut(i)
337	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
338	>	endif
339		endif
296	–	if (i_is_Sticky) then
297	–	thisRcut = getStickyCut(i)
298	–	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
299	–	endif
300	–	if (i_is_StickyP) then
301	–	thisRcut = getStickyPowerCut(i)
302	–	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
303	–	endif
304	–	if (i_is_GB) then
305	–	thisRcut = getGayBerneCut(i)
306	–	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
307	–	endif
308	–	if (i_is_EAM) then
309	–	thisRcut = getEAMCut(i)
310	–	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
311	–	endif
312	–	if (i_is_Shape) then
313	–	thisRcut = getShapeCut(i)
314	–	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
315	–	endif
340
341	+
342		if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
343		biggestAtypeCutoff = atypeMaxCutoff(i)
344		endif
345	+
346		endif
347		enddo
348
349	<	nGroupTypes = 0
349	>
350
351		istart = 1
352	+	jstart = 1
353		#ifdef IS_MPI
354		iend = nGroupsInRow
355	+	jend = nGroupsInCol
356		#else
357		iend = nGroups
358	+	jend = nGroups
359		#endif
360
361		!! allocate the groupToGtype and gtypeMaxCutoff here.
362	<	if(.not.allocated(groupToGtype)) then
363	<	allocate(groupToGtype(iend))
364	<	allocate(groupMaxCutoff(iend))
365	<	allocate(gtypeMaxCutoff(iend))
362	>	if(.not.allocated(groupToGtypeRow)) then
363	>	! allocate(groupToGtype(iend))
364	>	allocate(groupToGtypeRow(iend))
365	>	else
366	>	deallocate(groupToGtypeRow)
367	>	allocate(groupToGtypeRow(iend))
368		endif
369	+	if(.not.allocated(groupMaxCutoffRow)) then
370	+	allocate(groupMaxCutoffRow(iend))
371	+	else
372	+	deallocate(groupMaxCutoffRow)
373	+	allocate(groupMaxCutoffRow(iend))
374	+	end if
375	+
376	+	if(.not.allocated(gtypeMaxCutoffRow)) then
377	+	allocate(gtypeMaxCutoffRow(iend))
378	+	else
379	+	deallocate(gtypeMaxCutoffRow)
380	+	allocate(gtypeMaxCutoffRow(iend))
381	+	endif
382	+
383	+
384	+	#ifdef IS_MPI
385	+	! We only allocate new storage if we are in MPI because Ncol /= Nrow
386	+	if(.not.associated(groupToGtypeCol)) then
387	+	allocate(groupToGtypeCol(jend))
388	+	else
389	+	deallocate(groupToGtypeCol)
390	+	allocate(groupToGtypeCol(jend))
391	+	end if
392	+
393	+	if(.not.associated(groupToGtypeCol)) then
394	+	allocate(groupToGtypeCol(jend))
395	+	else
396	+	deallocate(groupToGtypeCol)
397	+	allocate(groupToGtypeCol(jend))
398	+	end if
399	+	if(.not.associated(gtypeMaxCutoffCol)) then
400	+	allocate(gtypeMaxCutoffCol(jend))
401	+	else
402	+	deallocate(gtypeMaxCutoffCol)
403	+	allocate(gtypeMaxCutoffCol(jend))
404	+	end if
405	+
406	+	groupMaxCutoffCol = 0.0_dp
407	+	gtypeMaxCutoffCol = 0.0_dp
408	+
409	+	#endif
410	+	groupMaxCutoffRow = 0.0_dp
411	+	gtypeMaxCutoffRow = 0.0_dp
412	+
413	+
414		!! first we do a single loop over the cutoff groups to find the
415		!! largest cutoff for any atypes present in this group. We also
416		!! create gtypes at this point.
417
418		tol = 1.0d-6
419	<
419	>	nGroupTypesRow = 0
420	>
421		do i = istart, iend
422		n_in_i = groupStartRow(i+1) - groupStartRow(i)
423	<	groupMaxCutoff(i) = 0.0_dp
423	>	groupMaxCutoffRow(i) = 0.0_dp
424		do ia = groupStartRow(i), groupStartRow(i+1)-1
425		atom1 = groupListRow(ia)
426		#ifdef IS_MPI
#	Line 351 \| Line 428 \| contains
428		#else
429		me_i = atid(atom1)
430		#endif
431	<	if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
432	<	groupMaxCutoff(i)=atypeMaxCutoff(me_i)
433	<	endif
431	>	if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
432	>	groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
433	>	endif
434		enddo
435	<	if (nGroupTypes.eq.0) then
436	<	nGroupTypes = nGroupTypes + 1
437	<	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
438	<	groupToGtype(i) = nGroupTypes
435	>
436	>	if (nGroupTypesRow.eq.0) then
437	>	nGroupTypesRow = nGroupTypesRow + 1
438	>	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
439	>	groupToGtypeRow(i) = nGroupTypesRow
440		else
441	<	do g = 1, nGroupTypes
442	<	if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).gt.tol) then
443	<	nGroupTypes = nGroupTypes + 1
444	<	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
445	<	groupToGtype(i) = nGroupTypes
368	<	else
369	<	groupToGtype(i) = g
441	>	GtypeFound = .false.
442	>	do g = 1, nGroupTypesRow
443	>	if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
444	>	groupToGtypeRow(i) = g
445	>	GtypeFound = .true.
446		endif
447		enddo
448	+	if (.not.GtypeFound) then
449	+	nGroupTypesRow = nGroupTypesRow + 1
450	+	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
451	+	groupToGtypeRow(i) = nGroupTypesRow
452	+	endif
453		endif
454	<	enddo
455	<
454	>	enddo
455	>
456	>	#ifdef IS_MPI
457	>	do j = jstart, jend
458	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
459	>	groupMaxCutoffCol(j) = 0.0_dp
460	>	do ja = groupStartCol(j), groupStartCol(j+1)-1
461	>	atom1 = groupListCol(ja)
462	>
463	>	me_j = atid_col(atom1)
464	>
465	>	if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
466	>	groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
467	>	endif
468	>	enddo
469	>
470	>	if (nGroupTypesCol.eq.0) then
471	>	nGroupTypesCol = nGroupTypesCol + 1
472	>	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
473	>	groupToGtypeCol(j) = nGroupTypesCol
474	>	else
475	>	GtypeFound = .false.
476	>	do g = 1, nGroupTypesCol
477	>	if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
478	>	groupToGtypeCol(j) = g
479	>	GtypeFound = .true.
480	>	endif
481	>	enddo
482	>	if (.not.GtypeFound) then
483	>	nGroupTypesCol = nGroupTypesCol + 1
484	>	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
485	>	groupToGtypeCol(j) = nGroupTypesCol
486	>	endif
487	>	endif
488	>	enddo
489	>
490	>	#else
491	>	! Set pointers to information we just found
492	>	nGroupTypesCol = nGroupTypesRow
493	>	groupToGtypeCol => groupToGtypeRow
494	>	gtypeMaxCutoffCol => gtypeMaxCutoffRow
495	>	groupMaxCutoffCol => groupMaxCutoffRow
496	>	#endif
497	>
498	>
499	>
500	>
501	>
502		!! allocate the gtypeCutoffMap here.
503	<	allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
503	>	allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
504		!! then we do a double loop over all the group TYPES to find the cutoff
505		!! map between groups of two types
506	<
507	<	do i = 1, nGroupTypes
508	<	do j = 1, nGroupTypes
506	>	tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
507	>
508	>	do i = 1, nGroupTypesRow
509	>	do j = 1, nGroupTypesCol
510
511		select case(cutoffPolicy)
512		case(TRADITIONAL_CUTOFF_POLICY)
513	<	thisRcut = maxval(gtypeMaxCutoff)
513	>	thisRcut = tradRcut
514		case(MIX_CUTOFF_POLICY)
515	<	thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
515	>	thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
516		case(MAX_CUTOFF_POLICY)
517	<	thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
517	>	thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
518		case default
519		call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
520		return
#	Line 394 \| Line 522 \| contains
522		gtypeCutoffMap(i,j)%rcut = thisRcut
523		gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
524		skin = defaultRlist - defaultRcut
525	+	listSkin = skin ! set neighbor list skin thickness
526		gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
527
528	+	! sanity check
529	+
530	+	if (haveDefaultCutoffs) then
531	+	if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
532	+	call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
533	+	endif
534	+	endif
535		enddo
536		enddo
537	+	if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
538	+	if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
539	+	if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
540	+	#ifdef IS_MPI
541	+	if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
542	+	if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
543	+	#endif
544	+	groupMaxCutoffCol => null()
545	+	gtypeMaxCutoffCol => null()
546
547		haveGtypeCutoffMap = .true.
548	<
404	<	end subroutine createGtypeCutoffMap
405	<
406	<	subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
407	<	real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
408	<	integer, intent(in) :: cutPolicy
409	<
410	<	defaultRcut = defRcut
411	<	defaultRsw = defRsw
412	<	defaultRlist = defRlist
413	<	cutoffPolicy = cutPolicy
414	<	end subroutine setDefaultCutoffs
415	<
416	<	subroutine setCutoffPolicy(cutPolicy)
548	>	end subroutine createGtypeCutoffMap
549
550	+	subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
551	+	real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
552		integer, intent(in) :: cutPolicy
553	+
554	+	defaultRcut = defRcut
555	+	defaultRsw = defRsw
556	+	defaultRlist = defRlist
557		cutoffPolicy = cutPolicy
420	–	call createGtypeCutoffMap()
558
559	+	haveDefaultCutoffs = .true.
560	+	end subroutine setDefaultCutoffs
561	+
562	+	subroutine setCutoffPolicy(cutPolicy)
563	+
564	+	integer, intent(in) :: cutPolicy
565	+	cutoffPolicy = cutPolicy
566	+	call createGtypeCutoffMap()
567		end subroutine setCutoffPolicy
568
569
570		subroutine setSimVariables()
571		SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
572		SIM_uses_EAM = SimUsesEAM()
428	–	SIM_uses_RF = SimUsesRF()
573		SIM_requires_postpair_calc = SimRequiresPostpairCalc()
574		SIM_requires_prepair_calc = SimRequiresPrepairCalc()
575		SIM_uses_PBC = SimUsesPBC()
#	Line 495 \| Line 639 \| contains
639		end subroutine doReadyCheck
640
641
642	<	subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
642	>	subroutine init_FF(thisESM, thisStat)
643
644	<	logical, intent(in) :: use_RF
501	<	logical, intent(in) :: use_UW
502	<	logical, intent(in) :: use_DW
644	>	integer, intent(in) :: thisESM
645		integer, intent(out) :: thisStat
646		integer :: my_status, nMatches
505	–	integer :: corrMethod
647		integer, pointer :: MatchList(:) => null()
648		real(kind=dp) :: rcut, rrf, rt, dielect
649
650		!! assume things are copacetic, unless they aren't
651		thisStat = 0
652
653	<	!! Fortran's version of a cast:
513	<	FF_uses_RF = use_RF
653	>	electrostaticSummationMethod = thisESM
654
515	–	!! set the electrostatic correction method
516	–	if (use_UW) then
517	–	corrMethod = 1
518	–	elseif (use_DW) then
519	–	corrMethod = 2
520	–	else
521	–	corrMethod = 0
522	–	endif
523	–
655		!! init_FF is called after all of the atom types have been
656		!! defined in atype_module using the new_atype subroutine.
657		!!
#	Line 550 \| Line 681 \| contains
681
682		haveSaneForceField = .true.
683
684	<	!! check to make sure the FF_uses_RF setting makes sense
684	>	!! check to make sure the reaction field setting makes sense
685
686		if (FF_uses_Dipoles) then
687	<	if (FF_uses_RF) then
687	>	if (electrostaticSummationMethod == REACTION_FIELD) then
688		dielect = getDielect()
689		call initialize_rf(dielect)
690		endif
691		else
692	<	if (FF_uses_RF) then
692	>	if (electrostaticSummationMethod == REACTION_FIELD) then
693		write(default_error,*) 'Using Reaction Field with no dipoles? Huh?'
694		thisStat = -1
695		haveSaneForceField = .false.
#	Line 618 \| Line 749 \| contains
749
750		!! Stress Tensor
751		real( kind = dp), dimension(9) :: tau
752	<	real ( kind = dp ) :: pot
752	>	real ( kind = dp ),dimension(POT_ARRAY_SIZE) :: pot
753		logical ( kind = 2) :: do_pot_c, do_stress_c
754		logical :: do_pot
755		logical :: do_stress
756		logical :: in_switching_region
757		#ifdef IS_MPI
758	<	real( kind = DP ) :: pot_local
758	>	real( kind = DP ), dimension(POT_ARRAY_SIZE) :: pot_local
759		integer :: nAtomsInRow
760		integer :: nAtomsInCol
761		integer :: nprocs
#	Line 649 \| Line 780 \| contains
780		integer :: propPack_i, propPack_j
781		integer :: loopStart, loopEnd, loop
782		integer :: iHash
783	<	real(kind=dp) :: listSkin = 1.0
783	>
784
785		!! initialize local variables
786
#	Line 774 \| Line 905 \| contains
905		me_j = atid(j)
906		call get_interatomic_vector(q_group(:,i), &
907		q_group(:,j), d_grp, rgrpsq)
908	<	#endif
908	>	#endif
909
910	<	if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
910	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
911		if (update_nlist) then
912		nlist = nlist + 1
913
#	Line 897 \| Line 1028 \| contains
1028		endif
1029		end if
1030		enddo
1031	+
1032		enddo outer
1033
1034		if (update_nlist) then
#	Line 956 \| Line 1088 \| contains
1088
1089		if (do_pot) then
1090		! scatter/gather pot_row into the members of my column
1091	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1092	<
1091	>	do i = 1,POT_ARRAY_SIZE
1092	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1093	>	end do
1094		! scatter/gather pot_local into all other procs
1095		! add resultant to get total pot
1096		do i = 1, nlocal
1097	<	pot_local = pot_local + pot_Temp(i)
1097	>	pot_local(1:POT_ARRAY_SIZE) = pot_local(1:POT_ARRAY_SIZE &
1098	>	+ pot_Temp(1:POT_ARRAY_SIZE,i)
1099		enddo
1100
1101		pot_Temp = 0.0_DP
1102	<
1103	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1102	>	do i = 1,POT_ARRAY_SIZE
1103	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1104	>	end do
1105		do i = 1, nlocal
1106	<	pot_local = pot_local + pot_Temp(i)
1106	>	pot_local(1:POT_ARRAY_SIZE) = pot_local(1:POT_ARRAY_SIZE)&
1107	>	+ pot_Temp(1:POT_ARRAY_SIZE,i)
1108		enddo
1109
1110		endif
#	Line 976 \| Line 1112 \| contains
1112
1113		if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1114
1115	<	if (FF_uses_RF .and. SIM_uses_RF) then
1115	>	if (electrostaticSummationMethod == REACTION_FIELD) then
1116
1117		#ifdef IS_MPI
1118		call scatter(rf_Row,rf,plan_atom_row_3d)
#	Line 1007 \| Line 1143 \| contains
1143		!! potential and torques:
1144		call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1145		#ifdef IS_MPI
1146	<	pot_local = pot_local + rfpot
1146	>	pot_local(RF_POT) = pot_local(RF_POT) + rfpot
1147		#else
1148	<	pot = pot + rfpot
1148	>	pot(RF_POT) = pot(RF_POT) + rfpot
1149
1150		#endif
1151		endif
#	Line 1021 \| Line 1157 \| contains
1157		#ifdef IS_MPI
1158
1159		if (do_pot) then
1160	<	pot = pot + pot_local
1160	>	pot(1:SIZE_POT_ARRAY) = pot(1:SIZE_POT_ARRAY) &
1161	>	+ pot_local(1:SIZE_POT_ARRAY)
1162		!! we assume the c code will do the allreduce to get the total potential
1163		!! we could do it right here if we needed to...
1164		endif
#	Line 1047 \| Line 1184 \| contains
1184		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1185		eFrame, A, f, t, pot, vpair, fpair)
1186
1187	<	real( kind = dp ) :: pot, vpair, sw
1187	>	real( kind = dp ) :: vpair, sw
1188	>	real( kind = dp ), dimension(POT_ARRAY_SIZE) :: pot
1189		real( kind = dp ), dimension(3) :: fpair
1190		real( kind = dp ), dimension(nLocal) :: mfact
1191		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1060 \| Line 1198 \| contains
1198		real ( kind = dp ), intent(inout) :: rijsq
1199		real ( kind = dp ) :: r
1200		real ( kind = dp ), intent(inout) :: d(3)
1063	–	real ( kind = dp ) :: ebalance
1201		integer :: me_i, me_j
1202
1203		integer :: iHash
#	Line 1080 \| Line 1217 \| contains
1217		iHash = InteractionHash(me_i, me_j)
1218
1219		if ( iand(iHash, LJ_PAIR).ne.0 ) then
1220	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1220	>	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot(LJ_POT), f, do_pot)
1221		endif
1222
1223		if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1224		call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1225	<	pot, eFrame, f, t, do_pot, corrMethod)
1225	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1226
1227	<	if (FF_uses_RF .and. SIM_uses_RF) then
1227	>	if (electrostaticSummationMethod == REACTION_FIELD) then
1228
1229		! CHECK ME (RF needs to know about all electrostatic types)
1230		call accumulate_rf(i, j, r, eFrame, sw)
#	Line 1098 \| Line 1235 \| contains
1235
1236		if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1237		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1238	<	pot, A, f, t, do_pot)
1238	>	pot(STICKY_POT), A, f, t, do_pot)
1239		endif
1240
1241		if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1242		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1243	<	pot, A, f, t, do_pot)
1243	>	pot(STICKYPOWER_POT), A, f, t, do_pot)
1244		endif
1245
1246		if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1247		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1248	<	pot, A, f, t, do_pot)
1248	>	pot(GAYBERNE_POT), A, f, t, do_pot)
1249		endif
1250
1251		if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1252		! call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1253	<	! pot, A, f, t, do_pot)
1253	>	! pot(GAYBERNE_LJ_POT), A, f, t, do_pot)
1254		endif
1255
1256		if ( iand(iHash, EAM_PAIR).ne.0 ) then
1257	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1257	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot(EAM_POT), f, &
1258		do_pot)
1259		endif
1260
1261		if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1262		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1263	<	pot, A, f, t, do_pot)
1263	>	pot(SHAPE_POT), A, f, t, do_pot)
1264		endif
1265
1266		if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1267		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1268	<	pot, A, f, t, do_pot)
1268	>	pot(SHAPE_LJ_POT), A, f, t, do_pot)
1269		endif
1270
1271		end subroutine do_pair
#	Line 1136 \| Line 1273 \| contains
1273		subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1274		do_pot, do_stress, eFrame, A, f, t, pot)
1275
1276	<	real( kind = dp ) :: pot, sw
1276	>	real( kind = dp ) :: sw
1277	>	real( kind = dp ), dimension(POT_ARRAY_SIZE) :: pot
1278		real( kind = dp ), dimension(9,nLocal) :: eFrame
1279		real (kind=dp), dimension(9,nLocal) :: A
1280		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1150 \| Line 1288 \| contains
1288
1289		integer :: me_i, me_j, iHash
1290
1291	+	r = sqrt(rijsq)
1292	+
1293		#ifdef IS_MPI
1294		me_i = atid_row(i)
1295		me_j = atid_col(j)
#	Line 1169 \| Line 1309 \| contains
1309
1310		subroutine do_preforce(nlocal,pot)
1311		integer :: nlocal
1312	<	real( kind = dp ) :: pot
1312	>	real( kind = dp ),dimension(POT_ARRAY_SIZE) :: pot
1313
1314		if (FF_uses_EAM .and. SIM_uses_EAM) then
1315	<	call calc_EAM_preforce_Frho(nlocal,pot)
1315	>	call calc_EAM_preforce_Frho(nlocal,pot(EAM_POT))
1316		endif
1317
1318
#	Line 1367 \| Line 1507 \| contains
1507
1508		function FF_RequiresPostpairCalc() result(doesit)
1509		logical :: doesit
1510	<	doesit = FF_uses_RF
1510	>	if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1511		end function FF_RequiresPostpairCalc
1512
1513		#ifdef PROFILE

Diff Legend

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

Comparing trunk/OOPSE-2.0/src/UseTheForce/doForces.F90 (file contents): Revision 2284 by gezelter, Wed Sep 7 19:18:54 2005 UTC vs. Revision 2356 by chuckv, Wed Oct 12 19:55:26 2005 UTC

Diff Legend

Comparing trunk/OOPSE-2.0/src/UseTheForce/doForces.F90 (file contents):
Revision 2284 by gezelter, Wed Sep 7 19:18:54 2005 UTC vs.
Revision 2356 by chuckv, Wed Oct 12 19:55:26 2005 UTC