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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2286 by gezelter, Wed Sep 7 22:08:39 2005 UTC vs.
Revision 2381 by chrisfen, Tue Oct 18 15:01:42 2005 UTC

#	Line 45 \| Line 45
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.39 2005-09-07 22:08:39 gezelter Exp $, $Date: 2005-09-07 22:08:39 $, $Name: not supported by cvs2svn $, $Revision: 1.39 $
48	>	!! @version $Id: doForces.F90,v 1.60 2005-10-18 15:01:42 chrisfen Exp $, $Date: 2005-10-18 15:01:42 $, $Name: not supported by cvs2svn $, $Revision: 1.60 $
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
62	<	use gb_pair
61	>	use gayberne
62		use shapes
63		use vector_class
64		use eam
#	Line 75 \| Line 74 \| module doForces
74		#include "UseTheForce/fSwitchingFunction.h"
75		#include "UseTheForce/fCutoffPolicy.h"
76		#include "UseTheForce/DarkSide/fInteractionMap.h"
77	+	#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
78
79
80		INTEGER, PARAMETER:: PREPAIR_LOOP = 1
#	Line 85 \| Line 85 \| module doForces
85		logical, save :: haveSaneForceField = .false.
86		logical, save :: haveInteractionHash = .false.
87		logical, save :: haveGtypeCutoffMap = .false.
88	+	logical, save :: haveDefaultCutoffs = .false.
89		logical, save :: haveRlist = .false.
90
91		logical, save :: FF_uses_DirectionalAtoms
92		logical, save :: FF_uses_Dipoles
93		logical, save :: FF_uses_GayBerne
94		logical, save :: FF_uses_EAM
94	–	logical, save :: FF_uses_RF
95
96		logical, save :: SIM_uses_DirectionalAtoms
97		logical, save :: SIM_uses_EAM
98	–	logical, save :: SIM_uses_RF
98		logical, save :: SIM_requires_postpair_calc
99		logical, save :: SIM_requires_prepair_calc
100		logical, save :: SIM_uses_PBC
101
102	<	integer, save :: corrMethod
102	>	integer, save :: electrostaticSummationMethod
103
104		public :: init_FF
105		public :: setDefaultCutoffs
#	Line 124 \| Line 123 \| module doForces
123		! Bit hash to determine pair-pair interactions.
124		integer, dimension(:,:), allocatable :: InteractionHash
125		real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
126	<	real(kind=dp), dimension(:), allocatable :: groupMaxCutoff
127	<	integer, dimension(:), allocatable :: groupToGtype
128	<	real(kind=dp), dimension(:), allocatable :: gtypeMaxCutoff
126	>	real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
127	>	real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
128	>
129	>	integer, dimension(:), allocatable, target :: groupToGtypeRow
130	>	integer, dimension(:), pointer :: groupToGtypeCol => null()
131	>
132	>	real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
133	>	real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
134		type ::gtypeCutoffs
135		real(kind=dp) :: rcut
136		real(kind=dp) :: rcutsq
#	Line 136 \| Line 140 \| module doForces
140
141		integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
142		real(kind=dp),save :: defaultRcut, defaultRsw, defaultRlist
143	+	real(kind=dp),save :: listSkin
144
145		contains
146
#	Line 179 \| Line 184 \| contains
184
185		if (.not. allocated(InteractionHash)) then
186		allocate(InteractionHash(nAtypes,nAtypes))
187	+	else
188	+	deallocate(InteractionHash)
189	+	allocate(InteractionHash(nAtypes,nAtypes))
190		endif
191
192		if (.not. allocated(atypeMaxCutoff)) then
193		allocate(atypeMaxCutoff(nAtypes))
194	+	else
195	+	deallocate(atypeMaxCutoff)
196	+	allocate(atypeMaxCutoff(nAtypes))
197		endif
198
199		do i = 1, nAtypes
#	Line 259 \| Line 270 \| contains
270		logical :: GtypeFound
271
272		integer :: myStatus, nAtypes, i, j, istart, iend, jstart, jend
273	<	integer :: n_in_i, me_i, ia, g, atom1, nGroupTypes
274	<	real(kind=dp):: thisSigma, bigSigma, thisRcut, tol, skin
273	>	integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
274	>	integer :: nGroupsInRow
275	>	integer :: nGroupsInCol
276	>	integer :: nGroupTypesRow,nGroupTypesCol
277	>	real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol, skin
278		real(kind=dp) :: biggestAtypeCutoff
279
280		stat = 0
#	Line 272 \| Line 286 \| contains
286		return
287		endif
288		endif
289	<
289	>	#ifdef IS_MPI
290	>	nGroupsInRow = getNgroupsInRow(plan_group_row)
291	>	nGroupsInCol = getNgroupsInCol(plan_group_col)
292	>	#endif
293		nAtypes = getSize(atypes)
294	<
294	>	! Set all of the initial cutoffs to zero.
295	>	atypeMaxCutoff = 0.0_dp
296		do i = 1, nAtypes
297		if (SimHasAtype(i)) then
298		call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
#	Line 285 \| Line 303 \| contains
303		call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
304		call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
305
306	<	atypeMaxCutoff(i) = 0.0_dp
307	<	if (i_is_LJ) then
308	<	thisRcut = getSigma(i) * 2.5_dp
309	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
310	<	endif
311	<	if (i_is_Elect) then
312	<	thisRcut = defaultRcut
313	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
314	<	endif
315	<	if (i_is_Sticky) then
316	<	thisRcut = getStickyCut(i)
317	<	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
306	>
307	>	if (haveDefaultCutoffs) then
308	>	atypeMaxCutoff(i) = defaultRcut
309	>	else
310	>	if (i_is_LJ) then
311	>	thisRcut = getSigma(i) * 2.5_dp
312	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
313	>	endif
314	>	if (i_is_Elect) then
315	>	thisRcut = defaultRcut
316	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
317	>	endif
318	>	if (i_is_Sticky) then
319	>	thisRcut = getStickyCut(i)
320	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
321	>	endif
322	>	if (i_is_StickyP) then
323	>	thisRcut = getStickyPowerCut(i)
324	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
325	>	endif
326	>	if (i_is_GB) then
327	>	thisRcut = getGayBerneCut(i)
328	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
329	>	endif
330	>	if (i_is_EAM) then
331	>	thisRcut = getEAMCut(i)
332	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
333	>	endif
334	>	if (i_is_Shape) then
335	>	thisRcut = getShapeCut(i)
336	>	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
337	>	endif
338		endif
301	–	if (i_is_StickyP) then
302	–	thisRcut = getStickyPowerCut(i)
303	–	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
304	–	endif
305	–	if (i_is_GB) then
306	–	thisRcut = getGayBerneCut(i)
307	–	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
308	–	endif
309	–	if (i_is_EAM) then
310	–	thisRcut = getEAMCut(i)
311	–	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
312	–	endif
313	–	if (i_is_Shape) then
314	–	thisRcut = getShapeCut(i)
315	–	if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
316	–	endif
339
340	+
341		if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
342		biggestAtypeCutoff = atypeMaxCutoff(i)
343		endif
344	+
345		endif
346		enddo
347
348	<	nGroupTypes = 0
348	>
349
350		istart = 1
351	+	jstart = 1
352		#ifdef IS_MPI
353		iend = nGroupsInRow
354	+	jend = nGroupsInCol
355		#else
356		iend = nGroups
357	+	jend = nGroups
358		#endif
359
360		!! allocate the groupToGtype and gtypeMaxCutoff here.
361	<	if(.not.allocated(groupToGtype)) then
362	<	allocate(groupToGtype(iend))
363	<	allocate(groupMaxCutoff(iend))
364	<	allocate(gtypeMaxCutoff(iend))
361	>	if(.not.allocated(groupToGtypeRow)) then
362	>	! allocate(groupToGtype(iend))
363	>	allocate(groupToGtypeRow(iend))
364	>	else
365	>	deallocate(groupToGtypeRow)
366	>	allocate(groupToGtypeRow(iend))
367		endif
368	+	if(.not.allocated(groupMaxCutoffRow)) then
369	+	allocate(groupMaxCutoffRow(iend))
370	+	else
371	+	deallocate(groupMaxCutoffRow)
372	+	allocate(groupMaxCutoffRow(iend))
373	+	end if
374	+
375	+	if(.not.allocated(gtypeMaxCutoffRow)) then
376	+	allocate(gtypeMaxCutoffRow(iend))
377	+	else
378	+	deallocate(gtypeMaxCutoffRow)
379	+	allocate(gtypeMaxCutoffRow(iend))
380	+	endif
381	+
382	+
383	+	#ifdef IS_MPI
384	+	! We only allocate new storage if we are in MPI because Ncol /= Nrow
385	+	if(.not.associated(groupToGtypeCol)) then
386	+	allocate(groupToGtypeCol(jend))
387	+	else
388	+	deallocate(groupToGtypeCol)
389	+	allocate(groupToGtypeCol(jend))
390	+	end if
391	+
392	+	if(.not.associated(groupToGtypeCol)) then
393	+	allocate(groupToGtypeCol(jend))
394	+	else
395	+	deallocate(groupToGtypeCol)
396	+	allocate(groupToGtypeCol(jend))
397	+	end if
398	+	if(.not.associated(gtypeMaxCutoffCol)) then
399	+	allocate(gtypeMaxCutoffCol(jend))
400	+	else
401	+	deallocate(gtypeMaxCutoffCol)
402	+	allocate(gtypeMaxCutoffCol(jend))
403	+	end if
404	+
405	+	groupMaxCutoffCol = 0.0_dp
406	+	gtypeMaxCutoffCol = 0.0_dp
407	+
408	+	#endif
409	+	groupMaxCutoffRow = 0.0_dp
410	+	gtypeMaxCutoffRow = 0.0_dp
411	+
412	+
413		!! first we do a single loop over the cutoff groups to find the
414		!! largest cutoff for any atypes present in this group. We also
415		!! create gtypes at this point.
416
417		tol = 1.0d-6
418	<
418	>	nGroupTypesRow = 0
419	>
420		do i = istart, iend
421		n_in_i = groupStartRow(i+1) - groupStartRow(i)
422	<	groupMaxCutoff(i) = 0.0_dp
422	>	groupMaxCutoffRow(i) = 0.0_dp
423		do ia = groupStartRow(i), groupStartRow(i+1)-1
424		atom1 = groupListRow(ia)
425		#ifdef IS_MPI
#	Line 352 \| Line 427 \| contains
427		#else
428		me_i = atid(atom1)
429		#endif
430	<	if (atypeMaxCutoff(me_i).gt.groupMaxCutoff(i)) then
431	<	groupMaxCutoff(i)=atypeMaxCutoff(me_i)
430	>	if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
431	>	groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
432		endif
433		enddo
434
435	<	if (nGroupTypes.eq.0) then
436	<	nGroupTypes = nGroupTypes + 1
437	<	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
438	<	groupToGtype(i) = nGroupTypes
435	>	if (nGroupTypesRow.eq.0) then
436	>	nGroupTypesRow = nGroupTypesRow + 1
437	>	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
438	>	groupToGtypeRow(i) = nGroupTypesRow
439		else
440		GtypeFound = .false.
441	<	do g = 1, nGroupTypes
442	<	if ( abs(groupMaxCutoff(i) - gtypeMaxCutoff(g)).lt.tol) then
443	<	groupToGtype(i) = g
441	>	do g = 1, nGroupTypesRow
442	>	if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
443	>	groupToGtypeRow(i) = g
444		GtypeFound = .true.
445		endif
446		enddo
447		if (.not.GtypeFound) then
448	<	nGroupTypes = nGroupTypes + 1
449	<	gtypeMaxCutoff(nGroupTypes) = groupMaxCutoff(i)
450	<	groupToGtype(i) = nGroupTypes
448	>	nGroupTypesRow = nGroupTypesRow + 1
449	>	gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
450	>	groupToGtypeRow(i) = nGroupTypesRow
451	>	endif
452	>	endif
453	>	enddo
454	>
455	>	#ifdef IS_MPI
456	>	do j = jstart, jend
457	>	n_in_j = groupStartCol(j+1) - groupStartCol(j)
458	>	groupMaxCutoffCol(j) = 0.0_dp
459	>	do ja = groupStartCol(j), groupStartCol(j+1)-1
460	>	atom1 = groupListCol(ja)
461	>
462	>	me_j = atid_col(atom1)
463	>
464	>	if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
465	>	groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
466	>	endif
467	>	enddo
468	>
469	>	if (nGroupTypesCol.eq.0) then
470	>	nGroupTypesCol = nGroupTypesCol + 1
471	>	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
472	>	groupToGtypeCol(j) = nGroupTypesCol
473	>	else
474	>	GtypeFound = .false.
475	>	do g = 1, nGroupTypesCol
476	>	if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
477	>	groupToGtypeCol(j) = g
478	>	GtypeFound = .true.
479	>	endif
480	>	enddo
481	>	if (.not.GtypeFound) then
482	>	nGroupTypesCol = nGroupTypesCol + 1
483	>	gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
484	>	groupToGtypeCol(j) = nGroupTypesCol
485		endif
486		endif
487		enddo
488
489	+	#else
490	+	! Set pointers to information we just found
491	+	nGroupTypesCol = nGroupTypesRow
492	+	groupToGtypeCol => groupToGtypeRow
493	+	gtypeMaxCutoffCol => gtypeMaxCutoffRow
494	+	groupMaxCutoffCol => groupMaxCutoffRow
495	+	#endif
496	+
497	+
498	+
499	+
500	+
501		!! allocate the gtypeCutoffMap here.
502	<	allocate(gtypeCutoffMap(nGroupTypes,nGroupTypes))
502	>	allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
503		!! then we do a double loop over all the group TYPES to find the cutoff
504		!! map between groups of two types
505	<
506	<	do i = 1, nGroupTypes
507	<	do j = 1, nGroupTypes
505	>	tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
506	>
507	>	do i = 1, nGroupTypesRow
508	>	do j = 1, nGroupTypesCol
509
510		select case(cutoffPolicy)
511		case(TRADITIONAL_CUTOFF_POLICY)
512	<	thisRcut = maxval(gtypeMaxCutoff)
512	>	thisRcut = tradRcut
513		case(MIX_CUTOFF_POLICY)
514	<	thisRcut = 0.5_dp * (gtypeMaxCutoff(i) + gtypeMaxCutoff(j))
514	>	thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
515		case(MAX_CUTOFF_POLICY)
516	<	thisRcut = max(gtypeMaxCutoff(i), gtypeMaxCutoff(j))
516	>	thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
517		case default
518		call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
519		return
#	Line 399 \| Line 521 \| contains
521		gtypeCutoffMap(i,j)%rcut = thisRcut
522		gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
523		skin = defaultRlist - defaultRcut
524	+	listSkin = skin ! set neighbor list skin thickness
525		gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skin)**2
526
527	+	! sanity check
528	+
529	+	if (haveDefaultCutoffs) then
530	+	if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
531	+	call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
532	+	endif
533	+	endif
534		enddo
535		enddo
536	+	if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
537	+	if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
538	+	if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
539	+	#ifdef IS_MPI
540	+	if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
541	+	if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
542	+	#endif
543	+	groupMaxCutoffCol => null()
544	+	gtypeMaxCutoffCol => null()
545
546		haveGtypeCutoffMap = .true.
547	<
409	<	end subroutine createGtypeCutoffMap
410	<
411	<	subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
412	<	real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
413	<	integer, intent(in) :: cutPolicy
414	<
415	<	defaultRcut = defRcut
416	<	defaultRsw = defRsw
417	<	defaultRlist = defRlist
418	<	cutoffPolicy = cutPolicy
419	<	end subroutine setDefaultCutoffs
420	<
421	<	subroutine setCutoffPolicy(cutPolicy)
547	>	end subroutine createGtypeCutoffMap
548
549	+	subroutine setDefaultCutoffs(defRcut, defRsw, defRlist, cutPolicy)
550	+	real(kind=dp),intent(in) :: defRcut, defRsw, defRlist
551		integer, intent(in) :: cutPolicy
552	+
553	+	defaultRcut = defRcut
554	+	defaultRsw = defRsw
555	+	defaultRlist = defRlist
556		cutoffPolicy = cutPolicy
425	–	call createGtypeCutoffMap()
557
558	+	haveDefaultCutoffs = .true.
559	+	end subroutine setDefaultCutoffs
560	+
561	+	subroutine setCutoffPolicy(cutPolicy)
562	+
563	+	integer, intent(in) :: cutPolicy
564	+	cutoffPolicy = cutPolicy
565	+	call createGtypeCutoffMap()
566		end subroutine setCutoffPolicy
567
568
569		subroutine setSimVariables()
570		SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
571		SIM_uses_EAM = SimUsesEAM()
433	–	SIM_uses_RF = SimUsesRF()
572		SIM_requires_postpair_calc = SimRequiresPostpairCalc()
573		SIM_requires_prepair_calc = SimRequiresPrepairCalc()
574		SIM_uses_PBC = SimUsesPBC()
#	Line 500 \| Line 638 \| contains
638		end subroutine doReadyCheck
639
640
641	<	subroutine init_FF(use_RF, use_UW, use_DW, thisStat)
641	>	subroutine init_FF(thisESM, thisStat)
642
643	<	logical, intent(in) :: use_RF
506	<	logical, intent(in) :: use_UW
507	<	logical, intent(in) :: use_DW
643	>	integer, intent(in) :: thisESM
644		integer, intent(out) :: thisStat
645		integer :: my_status, nMatches
510	–	integer :: corrMethod
646		integer, pointer :: MatchList(:) => null()
647		real(kind=dp) :: rcut, rrf, rt, dielect
648
649		!! assume things are copacetic, unless they aren't
650		thisStat = 0
651
652	<	!! Fortran's version of a cast:
518	<	FF_uses_RF = use_RF
652	>	electrostaticSummationMethod = thisESM
653
520	–	!! set the electrostatic correction method
521	–	if (use_UW) then
522	–	corrMethod = 1
523	–	elseif (use_DW) then
524	–	corrMethod = 2
525	–	else
526	–	corrMethod = 0
527	–	endif
528	–
654		!! init_FF is called after all of the atom types have been
655		!! defined in atype_module using the new_atype subroutine.
656		!!
#	Line 554 \| Line 679 \| contains
679
680
681		haveSaneForceField = .true.
557	–
558	–	!! check to make sure the FF_uses_RF setting makes sense
559	–
560	–	if (FF_uses_Dipoles) then
561	–	if (FF_uses_RF) then
562	–	dielect = getDielect()
563	–	call initialize_rf(dielect)
564	–	endif
565	–	else
566	–	if (FF_uses_RF) then
567	–	write(default_error,*) 'Using Reaction Field with no dipoles? Huh?'
568	–	thisStat = -1
569	–	haveSaneForceField = .false.
570	–	return
571	–	endif
572	–	endif
682
683		if (FF_uses_EAM) then
684		call init_EAM_FF(my_status)
#	Line 581 \| Line 690 \| contains
690		end if
691		endif
692
584	–	if (FF_uses_GayBerne) then
585	–	call check_gb_pair_FF(my_status)
586	–	if (my_status .ne. 0) then
587	–	thisStat = -1
588	–	haveSaneForceField = .false.
589	–	return
590	–	endif
591	–	endif
592	–
693		if (.not. haveNeighborList) then
694		!! Create neighbor lists
695		call expandNeighborList(nLocal, my_status)
#	Line 623 \| Line 723 \| contains
723
724		!! Stress Tensor
725		real( kind = dp), dimension(9) :: tau
726	<	real ( kind = dp ) :: pot
726	>	real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
727		logical ( kind = 2) :: do_pot_c, do_stress_c
728		logical :: do_pot
729		logical :: do_stress
730		logical :: in_switching_region
731		#ifdef IS_MPI
732	<	real( kind = DP ) :: pot_local
732	>	real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
733		integer :: nAtomsInRow
734		integer :: nAtomsInCol
735		integer :: nprocs
#	Line 654 \| Line 754 \| contains
754		integer :: propPack_i, propPack_j
755		integer :: loopStart, loopEnd, loop
756		integer :: iHash
757	<	real(kind=dp) :: listSkin = 1.0
757	>
758
759		!! initialize local variables
760
#	Line 779 \| Line 879 \| contains
879		me_j = atid(j)
880		call get_interatomic_vector(q_group(:,i), &
881		q_group(:,j), d_grp, rgrpsq)
882	<	#endif
882	>	#endif
883
884	<	if (rgrpsq < gtypeCutoffMap(groupToGtype(i),groupToGtype(j))%rListsq) then
884	>	if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
885		if (update_nlist) then
886		nlist = nlist + 1
887
#	Line 902 \| Line 1002 \| contains
1002		endif
1003		end if
1004		enddo
1005	+
1006		enddo outer
1007
1008		if (update_nlist) then
#	Line 961 \| Line 1062 \| contains
1062
1063		if (do_pot) then
1064		! scatter/gather pot_row into the members of my column
1065	<	call scatter(pot_Row, pot_Temp, plan_atom_row)
1066	<
1065	>	do i = 1,LR_POT_TYPES
1066	>	call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1067	>	end do
1068		! scatter/gather pot_local into all other procs
1069		! add resultant to get total pot
1070		do i = 1, nlocal
1071	<	pot_local = pot_local + pot_Temp(i)
1071	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1072	>	+ pot_Temp(1:LR_POT_TYPES,i)
1073		enddo
1074
1075		pot_Temp = 0.0_DP
1076	<
1077	<	call scatter(pot_Col, pot_Temp, plan_atom_col)
1076	>	do i = 1,LR_POT_TYPES
1077	>	call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1078	>	end do
1079		do i = 1, nlocal
1080	<	pot_local = pot_local + pot_Temp(i)
1080	>	pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1081	>	+ pot_Temp(1:LR_POT_TYPES,i)
1082		enddo
1083
1084		endif
#	Line 981 \| Line 1086 \| contains
1086
1087		if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1088
1089	<	if (FF_uses_RF .and. SIM_uses_RF) then
1089	>	if (electrostaticSummationMethod == REACTION_FIELD) then
1090
1091		#ifdef IS_MPI
1092		call scatter(rf_Row,rf,plan_atom_row_3d)
#	Line 1012 \| Line 1117 \| contains
1117		!! potential and torques:
1118		call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
1119		#ifdef IS_MPI
1120	<	pot_local = pot_local + rfpot
1120	>	pot_local(ELECTROSTATIC_POT) = pot_local(ELECTROSTATIC_POT) + rfpot
1121		#else
1122	<	pot = pot + rfpot
1122	>	pot(ELECTROSTATIC_POT) = pot(ELECTROSTATIC_POT) + rfpot
1123
1124		#endif
1125		endif
#	Line 1026 \| Line 1131 \| contains
1131		#ifdef IS_MPI
1132
1133		if (do_pot) then
1134	<	pot = pot + pot_local
1134	>	pot(1:LR_POT_TYPES) = pot(1:LR_POT_TYPES) &
1135	>	+ pot_local(1:LR_POT_TYPES)
1136		!! we assume the c code will do the allreduce to get the total potential
1137		!! we could do it right here if we needed to...
1138		endif
#	Line 1052 \| Line 1158 \| contains
1158		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1159		eFrame, A, f, t, pot, vpair, fpair)
1160
1161	<	real( kind = dp ) :: pot, vpair, sw
1161	>	real( kind = dp ) :: vpair, sw
1162	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1163		real( kind = dp ), dimension(3) :: fpair
1164		real( kind = dp ), dimension(nLocal) :: mfact
1165		real( kind = dp ), dimension(9,nLocal) :: eFrame
#	Line 1065 \| Line 1172 \| contains
1172		real ( kind = dp ), intent(inout) :: rijsq
1173		real ( kind = dp ) :: r
1174		real ( kind = dp ), intent(inout) :: d(3)
1068	–	real ( kind = dp ) :: ebalance
1175		integer :: me_i, me_j
1176
1177		integer :: iHash
#	Line 1085 \| Line 1191 \| contains
1191		iHash = InteractionHash(me_i, me_j)
1192
1193		if ( iand(iHash, LJ_PAIR).ne.0 ) then
1194	<	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
1194	>	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1195	>	pot(VDW_POT), f, do_pot)
1196		endif
1197
1198		if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1199		call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1200	<	pot, eFrame, f, t, do_pot, corrMethod)
1200	>	pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1201
1202	<	if (FF_uses_RF .and. SIM_uses_RF) then
1202	>	if (electrostaticSummationMethod == REACTION_FIELD) then
1203
1204		! CHECK ME (RF needs to know about all electrostatic types)
1205		call accumulate_rf(i, j, r, eFrame, sw)
#	Line 1103 \| Line 1210 \| contains
1210
1211		if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1212		call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1213	<	pot, A, f, t, do_pot)
1213	>	pot(HB_POT), A, f, t, do_pot)
1214		endif
1215
1216		if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1217		call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1218	<	pot, A, f, t, do_pot)
1218	>	pot(HB_POT), A, f, t, do_pot)
1219		endif
1220
1221		if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1222		call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1223	<	pot, A, f, t, do_pot)
1223	>	pot(VDW_POT), A, f, t, do_pot)
1224		endif
1225
1226		if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1227	<	! call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1228	<	! pot, A, f, t, do_pot)
1227	>	call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1228	>	pot(VDW_POT), A, f, t, do_pot)
1229		endif
1230
1231		if ( iand(iHash, EAM_PAIR).ne.0 ) then
1232	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
1233	<	do_pot)
1232	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1233	>	pot(METALLIC_POT), f, do_pot)
1234		endif
1235
1236		if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1237		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1238	<	pot, A, f, t, do_pot)
1238	>	pot(VDW_POT), A, f, t, do_pot)
1239		endif
1240
1241		if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1242		call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1243	<	pot, A, f, t, do_pot)
1243	>	pot(VDW_POT), A, f, t, do_pot)
1244		endif
1245
1246		end subroutine do_pair
#	Line 1141 \| Line 1248 \| contains
1248		subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1249		do_pot, do_stress, eFrame, A, f, t, pot)
1250
1251	<	real( kind = dp ) :: pot, sw
1251	>	real( kind = dp ) :: sw
1252	>	real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1253		real( kind = dp ), dimension(9,nLocal) :: eFrame
1254		real (kind=dp), dimension(9,nLocal) :: A
1255		real (kind=dp), dimension(3,nLocal) :: f
#	Line 1155 \| Line 1263 \| contains
1263
1264		integer :: me_i, me_j, iHash
1265
1266	+	r = sqrt(rijsq)
1267	+
1268		#ifdef IS_MPI
1269		me_i = atid_row(i)
1270		me_j = atid_col(j)
#	Line 1174 \| Line 1284 \| contains
1284
1285		subroutine do_preforce(nlocal,pot)
1286		integer :: nlocal
1287	<	real( kind = dp ) :: pot
1287	>	real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1288
1289		if (FF_uses_EAM .and. SIM_uses_EAM) then
1290	<	call calc_EAM_preforce_Frho(nlocal,pot)
1290	>	call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1291		endif
1292
1293
#	Line 1372 \| Line 1482 \| contains
1482
1483		function FF_RequiresPostpairCalc() result(doesit)
1484		logical :: doesit
1485	<	doesit = FF_uses_RF
1485	>	if (electrostaticSummationMethod == REACTION_FIELD) doesit = .true.
1486		end function FF_RequiresPostpairCalc
1487
1488		#ifdef PROFILE

Diff Legend

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

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents): Revision 2286 by gezelter, Wed Sep 7 22:08:39 2005 UTC vs. Revision 2381 by chrisfen, Tue Oct 18 15:01:42 2005 UTC

Diff Legend

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2286 by gezelter, Wed Sep 7 22:08:39 2005 UTC vs.
Revision 2381 by chrisfen, Tue Oct 18 15:01:42 2005 UTC