ViewVC Help

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

root/group/trunk/OOPSE/libmdtools/do_Forces.F90

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 597 by mmeineke, Mon Jul 14 21:28:54 2003 UTC vs.
Revision 1150 by gezelter, Fri May 7 21:35:05 2004 UTC

#	Line 4 | Line 4
4
5		!! @author Charles F. Vardeman II
6		!! @author Matthew Meineke
7	<	!! @version $Id: do_Forces.F90,v 1.18 2003-07-14 21:28:54 mmeineke Exp $, $Date: 2003-07-14 21:28:54 $, $Name: not supported by cvs2svn $, $Revision: 1.18 $
7	>	!! @version $Id: do_Forces.F90,v 1.54 2004-05-07 21:35:04 gezelter Exp $, $Date: 2004-05-07 21:35:04 $, $Name: not supported by cvs2svn $, $Revision: 1.54 $
8
9		module do_Forces
10		use force_globals
11		use simulation
12		use definitions
13		use atype_module
14	+	use switcheroo
15		use neighborLists
16		use lj
17		use sticky_pair
18		use dipole_dipole
19	+	use charge_charge
20		use reaction_field
21		use gb_pair
22	+	use vector_class
23	+	use eam
24	+	use status
25		#ifdef IS_MPI
26		use mpiSimulation
27		#endif
#	Line 26 | Line 31 | module do_Forces
31
32		#define __FORTRAN90
33		#include "fForceField.h"
34	+	#include "fSwitchingFunction.h"
35
36	<	logical, save :: do_forces_initialized = .false.
36	>	logical, save :: haveRlist = .false.
37	>	logical, save :: haveNeighborList = .false.
38	>	logical, save :: havePolicies = .false.
39	>	logical, save :: haveSIMvariables = .false.
40	>	logical, save :: havePropertyMap = .false.
41	>	logical, save :: haveSaneForceField = .false.
42		logical, save :: FF_uses_LJ
43		logical, save :: FF_uses_sticky
44	+	logical, save :: FF_uses_charges
45		logical, save :: FF_uses_dipoles
46		logical, save :: FF_uses_RF
47		logical, save :: FF_uses_GB
48		logical, save :: FF_uses_EAM
49	+	logical, save :: SIM_uses_LJ
50	+	logical, save :: SIM_uses_sticky
51	+	logical, save :: SIM_uses_charges
52	+	logical, save :: SIM_uses_dipoles
53	+	logical, save :: SIM_uses_RF
54	+	logical, save :: SIM_uses_GB
55	+	logical, save :: SIM_uses_EAM
56	+	logical, save :: SIM_requires_postpair_calc
57	+	logical, save :: SIM_requires_prepair_calc
58	+	logical, save :: SIM_uses_directional_atoms
59	+	logical, save :: SIM_uses_PBC
60	+	logical, save :: SIM_uses_molecular_cutoffs
61
62	+	real(kind=dp), save :: rlist, rlistsq
63	+
64		public :: init_FF
65		public :: do_force_loop
66	+	public :: setRlistDF
67
68	+	#ifdef PROFILE
69	+	public :: getforcetime
70	+	real, save :: forceTime = 0
71	+	real :: forceTimeInitial, forceTimeFinal
72	+	integer :: nLoops
73	+	#endif
74	+
75	+	type :: Properties
76	+	logical :: is_lj     = .false.
77	+	logical :: is_sticky = .false.
78	+	logical :: is_dp     = .false.
79	+	logical :: is_gb     = .false.
80	+	logical :: is_eam    = .false.
81	+	logical :: is_charge = .false.
82	+	real(kind=DP) :: charge = 0.0_DP
83	+	real(kind=DP) :: dipole_moment = 0.0_DP
84	+	end type Properties
85	+
86	+	type(Properties), dimension(:),allocatable :: PropertyMap
87	+
88		contains
89
90	+	subroutine setRlistDF( this_rlist )
91	+
92	+	real(kind=dp) :: this_rlist
93	+
94	+	rlist = this_rlist
95	+	rlistsq = rlist * rlist
96	+
97	+	haveRlist = .true.
98	+
99	+	end subroutine setRlistDF
100	+
101	+	subroutine createPropertyMap(status)
102	+	integer :: nAtypes
103	+	integer :: status
104	+	integer :: i
105	+	logical :: thisProperty
106	+	real (kind=DP) :: thisDPproperty
107	+
108	+	status = 0
109	+
110	+	nAtypes = getSize(atypes)
111	+
112	+	if (nAtypes == 0) then
113	+	status = -1
114	+	return
115	+	end if
116	+
117	+	if (.not. allocated(PropertyMap)) then
118	+	allocate(PropertyMap(nAtypes))
119	+	endif
120	+
121	+	do i = 1, nAtypes
122	+	call getElementProperty(atypes, i, "is_LJ", thisProperty)
123	+	PropertyMap(i)%is_LJ = thisProperty
124	+
125	+	call getElementProperty(atypes, i, "is_Charge", thisProperty)
126	+	PropertyMap(i)%is_Charge = thisProperty
127	+
128	+	if (thisProperty) then
129	+	call getElementProperty(atypes, i, "charge", thisDPproperty)
130	+	PropertyMap(i)%charge = thisDPproperty
131	+	endif
132	+
133	+	call getElementProperty(atypes, i, "is_DP", thisProperty)
134	+	PropertyMap(i)%is_DP = thisProperty
135	+
136	+	if (thisProperty) then
137	+	call getElementProperty(atypes, i, "dipole_moment", thisDPproperty)
138	+	PropertyMap(i)%dipole_moment = thisDPproperty
139	+	endif
140	+
141	+	call getElementProperty(atypes, i, "is_Sticky", thisProperty)
142	+	PropertyMap(i)%is_Sticky = thisProperty
143	+	call getElementProperty(atypes, i, "is_GB", thisProperty)
144	+	PropertyMap(i)%is_GB = thisProperty
145	+	call getElementProperty(atypes, i, "is_EAM", thisProperty)
146	+	PropertyMap(i)%is_EAM = thisProperty
147	+	end do
148	+
149	+	havePropertyMap = .true.
150	+
151	+	end subroutine createPropertyMap
152	+
153	+	subroutine setSimVariables()
154	+	SIM_uses_LJ = SimUsesLJ()
155	+	SIM_uses_sticky = SimUsesSticky()
156	+	SIM_uses_charges = SimUsesCharges()
157	+	SIM_uses_dipoles = SimUsesDipoles()
158	+	SIM_uses_RF = SimUsesRF()
159	+	SIM_uses_GB = SimUsesGB()
160	+	SIM_uses_EAM = SimUsesEAM()
161	+	SIM_requires_postpair_calc = SimRequiresPostpairCalc()
162	+	SIM_requires_prepair_calc = SimRequiresPrepairCalc()
163	+	SIM_uses_directional_atoms = SimUsesDirectionalAtoms()
164	+	SIM_uses_PBC = SimUsesPBC()
165	+	!SIM_uses_molecular_cutoffs = SimUsesMolecularCutoffs()
166	+
167	+	haveSIMvariables = .true.
168	+
169	+	return
170	+	end subroutine setSimVariables
171	+
172	+	subroutine doReadyCheck(error)
173	+	integer, intent(out) :: error
174	+
175	+	integer :: myStatus
176	+
177	+	error = 0
178	+
179	+	if (.not. havePropertyMap) then
180	+
181	+	myStatus = 0
182	+
183	+	call createPropertyMap(myStatus)
184	+
185	+	if (myStatus .ne. 0) then
186	+	write(default_error, *) 'createPropertyMap failed in do_Forces!'
187	+	error = -1
188	+	return
189	+	endif
190	+	endif
191	+
192	+	if (.not. haveSIMvariables) then
193	+	call setSimVariables()
194	+	endif
195	+
196	+	if (.not. haveRlist) then
197	+	write(default_error, *) 'rList has not been set in do_Forces!'
198	+	error = -1
199	+	return
200	+	endif
201	+
202	+	if (SIM_uses_LJ .and. FF_uses_LJ) then
203	+	if (.not. havePolicies) then
204	+	write(default_error, *) 'LJ mixing Policies have not been set in do_Forces!'
205	+	error = -1
206	+	return
207	+	endif
208	+	endif
209	+
210	+	if (.not. haveNeighborList) then
211	+	write(default_error, *) 'neighbor list has not been initialized in do_Forces!'
212	+	error = -1
213	+	return
214	+	end if
215	+
216	+	if (.not. haveSaneForceField) then
217	+	write(default_error, *) 'Force Field is not sane in do_Forces!'
218	+	error = -1
219	+	return
220	+	end if
221	+
222	+	#ifdef IS_MPI
223	+	if (.not. isMPISimSet()) then
224	+	write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
225	+	error = -1
226	+	return
227	+	endif
228	+	#endif
229	+	return
230	+	end subroutine doReadyCheck
231	+
232	+
233		subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
234
235		integer, intent(in) :: LJMIXPOLICY
#	Line 64 | Line 254 | contains
254
255		FF_uses_LJ = .false.
256		FF_uses_sticky = .false.
257	+	FF_uses_charges = .false.
258		FF_uses_dipoles = .false.
259		FF_uses_GB = .false.
260		FF_uses_EAM = .false.
261
262		call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList)
263		if (nMatches .gt. 0) FF_uses_LJ = .true.
264	<
264	>
265	>	call getMatchingElementList(atypes, "is_Charge", .true., nMatches, MatchList)
266	>	if (nMatches .gt. 0) FF_uses_charges = .true.
267	>
268		call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
269		if (nMatches .gt. 0) FF_uses_dipoles = .true.
270
#	Line 84 | Line 278 | contains
278		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
279		if (nMatches .gt. 0) FF_uses_EAM = .true.
280
281	+	!! Assume sanity (for the sake of argument)
282	+	haveSaneForceField = .true.
283	+	!!
284	+	if (FF_uses_charges) then
285	+	dielect = getDielect()
286	+	call initialize_charge(dielect)
287	+	endif
288	+
289	+
290		!! check to make sure the FF_uses_RF setting makes sense
291
292		if (FF_uses_dipoles) then
90	–	rrf = getRrf()
91	–	rt = getRt()
92	–	call initialize_dipole(rrf, rt)
293		if (FF_uses_RF) then
294		dielect = getDielect()
295	<	call initialize_rf(rrf, rt, dielect)
295	>	call initialize_rf(dielect)
296		endif
297		else
298		if (FF_uses_RF) then
299		write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
300		thisStat = -1
301	+	haveSaneForceField = .false.
302		return
303		endif
304	<	endif
304	>	endif
305
306		if (FF_uses_LJ) then
307
107	–	call getRcut(rcut)
108	–
308		select case (LJMIXPOLICY)
309		case (LB_MIXING_RULE)
310	<	call init_lj_FF(LB_MIXING_RULE, rcut, my_status)
310	>	call init_lj_FF(LB_MIXING_RULE, my_status)
311		case (EXPLICIT_MIXING_RULE)
312	<	call init_lj_FF(EXPLICIT_MIXING_RULE, rcut, my_status)
312	>	call init_lj_FF(EXPLICIT_MIXING_RULE, my_status)
313		case default
314		write(default_error,*) 'unknown LJ Mixing Policy!'
315		thisStat = -1
316	+	haveSaneForceField = .false.
317		return
318		end select
319		if (my_status /= 0) then
320		thisStat = -1
321	+	haveSaneForceField = .false.
322		return
323		end if
324	+	havePolicies = .true.
325		endif
326
327		if (FF_uses_sticky) then
328		call check_sticky_FF(my_status)
329		if (my_status /= 0) then
330		thisStat = -1
331	+	haveSaneForceField = .false.
332		return
333		end if
334		endif
335	<
335	>
336	>
337	>	if (FF_uses_EAM) then
338	>	call init_EAM_FF(my_status)
339	>	if (my_status /= 0) then
340	>	write(default_error, *) "init_EAM_FF returned a bad status"
341	>	thisStat = -1
342	>	haveSaneForceField = .false.
343	>	return
344	>	end if
345	>	endif
346	>
347		if (FF_uses_GB) then
348		call check_gb_pair_FF(my_status)
349		if (my_status .ne. 0) then
350		thisStat = -1
351	+	haveSaneForceField = .false.
352		return
353		endif
354		endif
355
356		if (FF_uses_GB .and. FF_uses_LJ) then
357		endif
358	<	if (.not. do_forces_initialized) then
358	>	if (.not. haveNeighborList) then
359		!! Create neighbor lists
360	<	call expandNeighborList(getNlocal(), my_status)
360	>	call expandNeighborList(nLocal, my_status)
361		if (my_Status /= 0) then
362		write(default_error,*) "SimSetup: ExpandNeighborList returned error."
363		thisStat = -1
364		return
365		endif
366	+	haveNeighborList = .true.
367		endif
368
369	<	do_forces_initialized = .true.
370	<
369	>
370	>
371		end subroutine init_FF
372
373
374		!! Does force loop over i,j pairs. Calls do_pair to calculates forces.
375		!------------------------------------------------------------->
376	<	subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
377	<	error)
376	>	subroutine do_force_loop(q, q_group, A, u_l, f, t, tau, pot, &
377	>	do_pot_c, do_stress_c, error)
378		!! Position array provided by C, dimensioned by getNlocal
379	<	real ( kind = dp ), dimension(3,getNlocal()) :: q
379	>	real ( kind = dp ), dimension(3,nLocal) :: q
380	>	!! molecular center-of-mass position array
381	>	real ( kind = dp ), dimension(3,nGroup) :: q_group
382		!! Rotation Matrix for each long range particle in simulation.
383	<	real( kind = dp), dimension(9,getNlocal()) :: A
383	>	real( kind = dp), dimension(9,nLocal) :: A
384		!! Unit vectors for dipoles (lab frame)
385	<	real( kind = dp ), dimension(3,getNlocal()) :: u_l
385	>	real( kind = dp ), dimension(3,nLocal) :: u_l
386		!! Force array provided by C, dimensioned by getNlocal
387	<	real ( kind = dp ), dimension(3,getNlocal()) :: f
387	>	real ( kind = dp ), dimension(3,nLocal) :: f
388		!! Torsion array provided by C, dimensioned by getNlocal
389	<	real( kind = dp ), dimension(3,getNlocal()) :: t
389	>	real( kind = dp ), dimension(3,nLocal) :: t
390	>
391		!! Stress Tensor
392		real( kind = dp), dimension(9) :: tau
393		real ( kind = dp ) :: pot
394		logical ( kind = 2) :: do_pot_c, do_stress_c
395		logical :: do_pot
396		logical :: do_stress
397	+	logical :: in_switching_region
398		#ifdef IS_MPI
399		real( kind = DP ) :: pot_local
400		integer :: nrow
401		integer :: ncol
402	+	integer :: nprocs
403	+	integer :: nrow_group
404	+	integer :: ncol_group
405		#endif
183	–	integer :: nlocal
406		integer :: natoms
407		logical :: update_nlist
408		integer :: i, j, jbeg, jend, jnab
409	+	integer :: ia, jb, atom1, atom2
410		integer :: nlist
411	<	real( kind = DP ) ::  rijsq, rlistsq, rcutsq, rlist, rcut
412	<	real(kind=dp),dimension(3) :: d
411	>	real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vab
412	>	real( kind = DP ) :: sw, dswdr, swderiv, mf
413	>	real(kind=dp),dimension(3) :: d_atm, d_grp
414		real(kind=dp) :: rfpot, mu_i, virial
415	<	integer :: me_i
415	>	integer :: me_i, me_j
416		logical :: is_dp_i
417		integer :: neighborListSize
418		integer :: listerror, error
419		integer :: localError
420	+	integer :: propPack_i, propPack_j
421
422	+	real(kind=dp) :: listSkin = 1.0
423	+
424		!! initialize local variables
425	<
425	>
426		#ifdef IS_MPI
427		pot_local = 0.0_dp
201	–	nlocal = getNlocal()
428		nrow   = getNrow(plan_row)
429		ncol   = getNcol(plan_col)
430	+	nrow_group   = getNrowGroup(plan_row)
431	+	ncol_group   = getNcolGroup(plan_col)
432		#else
205	–	nlocal = getNlocal()
433		natoms = nlocal
434		#endif
208	–
209	–	call getRcut(rcut,rc2=rcutsq)
210	–	call getRlist(rlist,rlistsq)
435
436	<	call check_initialization(localError)
436	>	call doReadyCheck(localError)
437		if ( localError .ne. 0 ) then
438	+	call handleError("do_force_loop", "Not Initialized")
439		error = -1
440		return
441		end if
442		call zero_work_arrays()
443	<
443	>
444		do_pot = do_pot_c
445		do_stress = do_stress_c
446	<
446	>
447		! Gather all information needed by all force loops:
448
449		#ifdef IS_MPI
450	+
451	+	call gather(q, q_Row, plan_row3d)
452	+	call gather(q, q_Col, plan_col3d)
453
454	<	call gather(q,q_Row,plan_row3d)
455	<	call gather(q,q_Col,plan_col3d)
454	>	call gather(q_group, q_group_Row, plan_row_Group_3d)
455	>	call gather(q_group, q_group_Col, plan_col_Group_3d)
456
457	<	if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
457	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
458		call gather(u_l,u_l_Row,plan_row3d)
459		call gather(u_l,u_l_Col,plan_col3d)
460
#	Line 236 | Line 464 | contains
464
465		#endif
466
467	<	if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
467	>	!! Begin force loop timing:
468	>	#ifdef PROFILE
469	>	call cpu_time(forceTimeInitial)
470	>	nloops = nloops + 1
471	>	#endif
472	>
473	>	if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
474		!! See if we need to update neighbor lists
475	<	call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)
475	>
476	>	call checkNeighborList(nGroup, q_group, listSkin, update_nlist)
477	>
478		!! if_mpi_gather_stuff_for_prepair
479		!! do_prepair_loop_if_needed
480		!! if_mpi_scatter_stuff_from_prepair
481		!! if_mpi_gather_stuff_from_prepair_to_main_loop
482	<	else
483	<	!! See if we need to update neighbor lists
484	<	call checkNeighborList(nlocal, q, rcut, rlist, update_nlist)
482	>
483	>	!--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
484	>	#ifdef IS_MPI
485	>
486	>	if (update_nlist) then
487	>
488	>	!! save current configuration, construct neighbor list,
489	>	!! and calculate forces
490	>
491	>	call saveNeighborList(nGroup, q_group)
492	>
493	>	neighborListSize = size(list)
494	>	nlist = 0
495	>
496	>	do i = 1, nrow_group
497	>	point(i) = nlist + 1
498	>
499	>	do j = 1, ncol_group
500	>
501	>	call get_interatomic_vector(q_group_Row(:,i), &
502	>	q_group_Col(:,j), d_grp, rgrpsq)
503	>
504	>	if (rgrpsq < rlistsq) then
505	>	nlist = nlist + 1
506	>
507	>	if (nlist > neighborListSize) then
508	>	call expandNeighborList(nGroup, listerror)
509	>	if (listerror /= 0) then
510	>	error = -1
511	>	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
512	>	return
513	>	end if
514	>	neighborListSize = size(list)
515	>	endif
516	>
517	>	list(nlist) = j
518	>
519	>	do ia = groupStart(i), groupStart(i+1)-1
520	>	atom1 = groupList(ia)
521	>
522	>	prepair_inner: do jb = groupStart(j), groupStart(j+1)-1
523	>	atom2 = groupList(jb)
524	>
525	>	if (skipThisPair(atom1, atom2)) cycle prepair_inner
526	>
527	>	call get_interatomic_vector(q_Row(:,atom1), &
528	>	q_Col(:,atom2), d_atm, ratmsq)
529	>
530	>	call do_prepair(atom1, atom2, ratmsq, d_atm, &
531	>	rgrpsq, d_grp, do_pot, do_stress, &
532	>	u_l, A, f, t, pot_local)
533	>
534	>	enddo prepair_inner
535	>	enddo
536	>	end if
537	>	enddo
538	>	enddo
539	>	point(nrow_group + 1) = nlist + 1
540	>
541	>	else  !! (of update_check)
542	>
543	>	! use the list to find the neighbors
544	>	do i = 1, nrow_group
545	>	JBEG = POINT(i)
546	>	JEND = POINT(i+1) - 1
547	>	! check that group i has neighbors
548	>	if (jbeg .le. jend) then
549	>
550	>	do jnab = jbeg, jend
551	>	j = list(jnab)
552	>
553	>	do ia = groupStart(i), groupStart(i+1)-1
554	>	atom1 = groupList(ia)
555	>
556	>	do jb = groupStart(j), groupStart(j+1)-1
557	>	atom2 = groupList(jb)
558	>
559	>	call get_interatomic_vector(q_Row(:,atom1), &
560	>	q_Col(:,atom2), d_atm, ratmsq)
561	>
562	>	call do_prepair(atom1, atom2, ratmsq, d_atm, &
563	>	rgrpsq, d_grp, do_pot, do_stress, &
564	>	u_l, A, f, t, pot_local)
565	>
566	>	enddo
567	>	enddo
568	>	enddo
569	>	endif
570	>	enddo
571	>	endif
572	>
573	>	#else
574	>
575	>	if (update_nlist) then
576	>
577	>	!! save current configuration, construct neighbor list,
578	>	!! and calculate forces
579	>
580	>	call saveNeighborList(nGroup, q_group)
581	>
582	>	neighborListSize = size(list)
583	>	nlist = 0
584	>
585	>	do i = 1, nGroup-1
586	>	point(i) = nlist + 1
587	>
588	>	do j = i+1, nGroup
589	>
590	>	call get_interatomic_vector(q_group(:,i), &
591	>	q_group(:,j), d_grp, rgrpsq)
592	>
593	>	if (rgrpsq < rlistsq) then
594	>	nlist = nlist + 1
595	>
596	>	if (nlist > neighborListSize) then
597	>	call expandNeighborList(nGroup, listerror)
598	>	if (listerror /= 0) then
599	>	error = -1
600	>	write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
601	>	return
602	>	end if
603	>	neighborListSize = size(list)
604	>	endif
605	>
606	>	list(nlist) = j
607	>
608	>	do ia = groupStart(i), groupStart(i+1)-1
609	>	atom1 = groupList(ia)
610	>
611	>	prepair_inner: do jb = groupStart(j), groupStart(j+1)-1
612	>	atom2 = groupList(jb)
613	>
614	>	if (skipThisPair(atom1, atom2)) cycle prepair_inner
615	>
616	>	call get_interatomic_vector(q(:,atom1), &
617	>	q(:,atom2), d_atm, ratmsq)
618	>
619	>	call do_prepair(atom1, atom2, ratmsq, d_atm, &
620	>	rgrpsq, d_grp, do_pot, do_stress, &
621	>	u_l, A, f, t, pot)
622	>
623	>	enddo prepair_inner
624	>	enddo
625	>	end if
626	>	enddo
627	>	enddo
628	>	point(nGroup) = nlist + 1
629	>
630	>	else  !! (of update_check)
631	>
632	>	! use the list to find the neighbors
633	>	do i = 1, nGroup-1
634	>	JBEG = POINT(i)
635	>	JEND = POINT(i+1) - 1
636	>	! check that group i has neighbors
637	>	if (jbeg .le. jend) then
638	>
639	>	do jnab = jbeg, jend
640	>	j = list(jnab)
641	>
642	>	do ia = groupStart(i), groupStart(i+1)-1
643	>	atom1 = groupList(ia)
644	>
645	>	do jb = groupStart(j), groupStart(j+1)-1
646	>	atom2 = groupList(jb)
647	>
648	>	call get_interatomic_vector(q(:,atom1), &
649	>	q(:,atom2), d_atm, ratmsq)
650	>
651	>	call do_prepair(atom1, atom2, ratmsq, d_atm, &
652	>	rgrpsq, d_grp, do_pot, do_stress, &
653	>	u_l, A, f, t, pot)
654	>
655	>	enddo
656	>	enddo
657	>	enddo
658	>	endif
659	>	enddo
660	>	endif
661	>
662	>	#endif
663	>
664	>	!! Do rest of preforce calculations
665	>	!! do necessary preforce calculations
666	>	call do_preforce(nlocal,pot)
667	>	! we have already updated the neighbor list set it to false...
668	>	update_nlist = .false.
669	>	else
670	>	!! See if we need to update neighbor lists for non pre-pair
671	>	call checkNeighborList(nGroup, q_group, listSkin, update_nlist)
672		endif
673
674	+	!---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>
675	+
676		#ifdef IS_MPI
677
678		if (update_nlist) then
679
680		!! save current configuration, construct neighbor list,
681		!! and calculate forces
257	–	call saveNeighborList(nlocal, q)
682
683	+	call saveNeighborList(nGroup, q_group)
684	+
685		neighborListSize = size(list)
686		nlist = 0
687
688	<	do i = 1, nrow
688	>	do i = 1, nrow_group
689		point(i) = nlist + 1
690
691	<	inner: do j = 1, ncol
691	>	do j = 1, ncol_group
692
693	<	if (skipThisPair(i,j)) cycle inner
693	>	call get_interatomic_vector(q_group_Row(:,i), &
694	>	q_group_Col(:,j), d_grp, rgrpsq)
695
696	<	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
270	<
271	<	if (rijsq <  rlistsq) then
272	<
696	>	if (rgrpsq < rlistsq) then
697		nlist = nlist + 1
698
699		if (nlist > neighborListSize) then
700	<	call expandNeighborList(nlocal, listerror)
700	>	call expandNeighborList(nGroup, listerror)
701		if (listerror /= 0) then
702		error = -1
703		write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
#	Line 283 | Line 707 | contains
707		endif
708
709		list(nlist) = j
710	<
711	<	if (rijsq <  rcutsq) then
712	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
713	<	u_l, A, f, t, pot_local)
710	>
711	>	vab = 0.0d0
712	>	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
713	>	in_switching_region)
714	>
715	>	do ia = groupStart(i), groupStart(i+1)-1
716	>	atom1 = groupList(ia)
717	>
718	>	inner: do jb = groupStart(j), groupStart(j+1)-1
719	>	atom2 = groupList(jb)
720	>
721	>	if (skipThisPair(atom1, atom2)) cycle inner
722	>
723	>	call get_interatomic_vector(q_Row(:,atom1), &
724	>	q_Col(:,atom2), d_atm, ratmsq)
725	>
726	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
727	>	do_pot, do_stress, &
728	>	u_l, A, f, t, pot_local, vpair)
729	>
730	>	vab = vab + vpair
731	>	enddo inner
732	>	enddo
733	>
734	>	if (in_switching_region) then
735	>	swderiv = vab*dswdr/rgrp
736	>
737	>	do ia=groupStart(i), groupStart(i+1)-1
738	>	atom1=groupList(ia)
739	>	mf = mfact(atom1)
740	>	f_Row(1,atom1) = f_Row(1,atom1) - swderiv*d_grp(1)*mf
741	>	f_Row(2,atom1) = f_Row(2,atom1) - swderiv*d_grp(2)*mf
742	>	f_Row(3,atom1) = f_Row(3,atom1) - swderiv*d_grp(3)*mf
743	>	enddo
744	>
745	>	do jb=groupStart(j), groupStart(j+1)-1
746	>	atom2=groupList(jb)
747	>	mf = mfact(atom2)
748	>	f_Col(1,atom2) = f_Col(1,atom2) + swderiv*d_grp(1)*mf
749	>	f_Col(2,atom2) = f_Col(2,atom2) + swderiv*d_grp(2)*mf
750	>	f_Col(3,atom2) = f_Col(3,atom2) + swderiv*d_grp(3)*mf
751	>	enddo
752		endif
291	–	endif
292	–	enddo inner
293	–	enddo
753
754	<	point(nrow + 1) = nlist + 1
754	>	end if
755	>	enddo
756	>	enddo
757	>	point(nrow_group + 1) = nlist + 1
758
759		else  !! (of update_check)
760	<
760	>
761		! use the list to find the neighbors
762	<	do i = 1, nrow
762	>	do i = 1, nrow_group
763		JBEG = POINT(i)
764		JEND = POINT(i+1) - 1
765	<	! check thiat molecule i has neighbors
765	>	! check that group i has neighbors
766		if (jbeg .le. jend) then
767
768		do jnab = jbeg, jend
769		j = list(jnab)
770	+
771	+	call get_interatomic_vector(q_group_Row(:,i), &
772	+	q_group_Col(:,j), d_grp, rgrpsq)
773	+
774	+	vab = 0.0d0
775	+	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
776	+	in_switching_region)
777	+
778	+	do ia = groupStart(i), groupStart(i+1)-1
779	+	atom1 = groupList(ia)
780	+
781	+	do jb = groupStart(j), groupStart(j+1)-1
782	+	atom2 = groupList(jb)
783	+
784	+	call get_interatomic_vector(q_Row(:,atom1), &
785	+	q_Col(:,atom2), d_atm, ratmsq)
786	+
787	+	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
788	+	do_pot, do_stress, &
789	+	u_l, A, f, t, pot_local, vpair)
790	+
791	+	vab = vab + vpair
792	+
793	+	enddo
794	+	enddo
795	+
796	+	if (in_switching_region) then
797	+	swderiv = vab*dswdr/rgrp
798	+
799	+	do ia=groupStart(i), groupStart(i+1)-1
800	+	atom1=groupList(ia)
801	+	mf = mfact(atom1)
802	+	f_Row(1,atom1) = f_Row(1,atom1) - swderiv*d_grp(1)*mf
803	+	f_Row(2,atom1) = f_Row(2,atom1) - swderiv*d_grp(2)*mf
804	+	f_Row(3,atom1) = f_Row(3,atom1) - swderiv*d_grp(3)*mf
805	+	enddo
806	+
807	+	do jb=groupStart(j), groupStart(j+1)-1
808	+	atom2=groupList(jb)
809	+	mf = mfact(atom2)
810	+	f_Col(1,atom2) = f_Col(1,atom2) + swderiv*d_grp(1)*mf
811	+	f_Col(2,atom2) = f_Col(2,atom2) + swderiv*d_grp(2)*mf
812	+	f_Col(3,atom2) = f_Col(3,atom2) + swderiv*d_grp(3)*mf
813	+	enddo
814	+	endif
815
309	–	call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
310	–	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
311	–	u_l, A, f, t, pot_local)
312	–
816		enddo
817		endif
818		enddo
819		endif
820
821		#else
822	<
822	>
823		if (update_nlist) then
824
825	<	! save current configuration, contruct neighbor list,
826	<	! and calculate forces
324	<	call saveNeighborList(natoms, q)
825	>	!! save current configuration, construct neighbor list,
826	>	!! and calculate forces
827
828	+	call saveNeighborList(nGroup, q_group)
829	+
830		neighborListSize = size(list)
831	<
328	<	nlist = 0
831	>	nlist = 0
832
833	<	do i = 1, natoms-1
833	>	do i = 1, nGroup-1
834		point(i) = nlist + 1
835
836	<	inner: do j = i+1, natoms
836	>	do j = i+1, nGroup
837
838	<	if (skipThisPair(i,j))  cycle inner
839	<
840	<	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
841	<
339	<
340	<	if (rijsq <  rlistsq) then
341	<
838	>	call get_interatomic_vector(q_group(:,i), &
839	>	q_group(:,j), d_grp, rgrpsq)
840	>
841	>	if (rgrpsq < rlistsq) then
842		nlist = nlist + 1
843	<
843	>
844		if (nlist > neighborListSize) then
845	<	call expandNeighborList(natoms, listerror)
845	>	call expandNeighborList(nGroup, listerror)
846		if (listerror /= 0) then
847		error = -1
848		write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
#	Line 352 | Line 852 | contains
852		endif
853
854		list(nlist) = j
855	+
856	+	vab = 0.0d0
857	+	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
858	+	in_switching_region)
859
860	<	if (rijsq <  rcutsq) then
861	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
862	<	u_l, A, f, t, pot)
860	>	do ia = groupStart(i), groupStart(i+1)-1
861	>	atom1 = groupList(ia)
862	>
863	>	inner: do jb = groupStart(j), groupStart(j+1)-1
864	>	atom2 = groupList(jb)
865	>
866	>	if (skipThisPair(atom1, atom2)) cycle inner
867	>
868	>	call get_interatomic_vector(q(:,atom1), &
869	>	q(:,atom2), d_atm, ratmsq)
870	>
871	>	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
872	>	do_pot, do_stress, &
873	>	u_l, A, f, t, pot, vpair)
874	>
875	>	vab = vab + vpair
876	>
877	>	enddo inner
878	>	enddo
879	>
880	>	if (in_switching_region) then
881	>	swderiv = vab*dswdr/rgrp
882	>	do ia=groupStart(i), groupStart(i+1)-1
883	>	atom1=groupList(ia)
884	>	mf = mfact(atom1)
885	>	f(1,atom1) = f(1,atom1) - swderiv*d_grp(1)*mf
886	>	f(2,atom1) = f(2,atom1) - swderiv*d_grp(2)*mf
887	>	f(3,atom1) = f(3,atom1) - swderiv*d_grp(3)*mf
888	>	enddo
889	>
890	>	do jb=groupStart(j), groupStart(j+1)-1
891	>	atom2=groupList(jb)
892	>	mf = mfact(atom2)
893	>	f(1,atom2) = f(1,atom2) + swderiv*d_grp(1)*mf
894	>	f(2,atom2) = f(2,atom2) + swderiv*d_grp(2)*mf
895	>	f(3,atom2) = f(3,atom2) + swderiv*d_grp(3)*mf
896	>	enddo
897		endif
898	<	endif
899	<	enddo inner
898	>
899	>	end if
900	>	enddo
901		enddo
902	+	point(nGroup) = nlist + 1
903
904	<	point(natoms) = nlist + 1
904	>	else  !! (of update_check)
905
366	–	else !! (update)
367	–
906		! use the list to find the neighbors
907	<	do i = 1, natoms-1
907	>	do i = 1, nGroup-1
908		JBEG = POINT(i)
909		JEND = POINT(i+1) - 1
910	<	! check thiat molecule i has neighbors
910	>	! check that group i has neighbors
911		if (jbeg .le. jend) then
912
913		do jnab = jbeg, jend
914		j = list(jnab)
915
916	<	call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
917	<	call do_pair(i, j, rijsq, d, do_pot, do_stress, &
380	<	u_l, A, f, t, pot)
916	>	call get_interatomic_vector(q_group(:,i), &
917	>	q_group(:,j), d_grp, rgrpsq)
918
919	+	vab = 0.0d0
920	+	call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
921	+	in_switching_region)
922	+
923	+	do ia = groupStart(i), groupStart(i+1)-1
924	+	atom1 = groupList(ia)
925	+
926	+	do jb = groupStart(j), groupStart(j+1)-1
927	+	atom2 = groupList(jb)
928	+
929	+	call get_interatomic_vector(q(:,atom1), &
930	+	q(:,atom2), d_atm, ratmsq)
931	+
932	+	call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
933	+	do_pot, do_stress, &
934	+	u_l, A, f, t, pot, vpair)
935	+
936	+	vab = vab + vpair
937	+
938	+	enddo
939	+	enddo
940	+
941	+	if (in_switching_region) then
942	+	swderiv = vab*dswdr/rgrp
943	+
944	+	do ia=groupStart(i), groupStart(i+1)-1
945	+	atom1=groupList(ia)
946	+	mf = mfact(atom1)
947	+	f(1,atom1) = f(1,atom1) - swderiv*d_grp(1)*mf
948	+	f(2,atom1) = f(2,atom1) - swderiv*d_grp(2)*mf
949	+	f(3,atom1) = f(3,atom1) - swderiv*d_grp(3)*mf
950	+	enddo
951	+
952	+	do jb=groupStart(j), groupStart(j+1)-1
953	+	atom2=groupList(jb)
954	+	mf = mfact(atom2)
955	+	f(1,atom2) = f(1,atom2) + swderiv*d_grp(1)*mf
956	+	f(2,atom2) = f(2,atom2) + swderiv*d_grp(2)*mf
957	+	f(3,atom2) = f(3,atom2) + swderiv*d_grp(3)*mf
958	+	enddo
959	+	endif
960		enddo
961		endif
962		enddo
963		endif
964	<
964	>
965		#endif
966
967		! phew, done with main loop.
968
969	+	!! Do timing
970	+	#ifdef PROFILE
971	+	call cpu_time(forceTimeFinal)
972	+	forceTime = forceTime + forceTimeFinal - forceTimeInitial
973	+	#endif
974	+
975		#ifdef IS_MPI
976		!!distribute forces
977	<
977	>
978		f_temp = 0.0_dp
979		call scatter(f_Row,f_temp,plan_row3d)
980		do i = 1,nlocal
981		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
982		end do
983	<
983	>
984		f_temp = 0.0_dp
985		call scatter(f_Col,f_temp,plan_col3d)
986		do i = 1,nlocal
987		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
988		end do
989
990	<	if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
990	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
991		t_temp = 0.0_dp
992		call scatter(t_Row,t_temp,plan_row3d)
993		do i = 1,nlocal
#	Line 420 | Line 1004 | contains
1004		if (do_pot) then
1005		! scatter/gather pot_row into the members of my column
1006		call scatter(pot_Row, pot_Temp, plan_row)
1007	<
1007	>
1008		! scatter/gather pot_local into all other procs
1009		! add resultant to get total pot
1010		do i = 1, nlocal
#	Line 428 | Line 1012 | contains
1012		enddo
1013
1014		pot_Temp = 0.0_DP
1015	<
1015	>
1016		call scatter(pot_Col, pot_Temp, plan_col)
1017		do i = 1, nlocal
1018		pot_local = pot_local + pot_Temp(i)
1019		enddo
1020	<
1021	<	endif
1020	>
1021	>	endif
1022		#endif
1023	<
1024	<	if (FF_RequiresPostpairCalc() .and. SimRequiresPostpairCalc()) then
1023	>
1024	>	if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
1025
1026	<	if (FF_uses_RF .and. SimUsesRF()) then
1026	>	if (FF_uses_RF .and. SIM_uses_RF) then
1027
1028		#ifdef IS_MPI
1029		call scatter(rf_Row,rf,plan_row3d)
#	Line 449 | Line 1033 | contains
1033		end do
1034		#endif
1035
1036	<	do i = 1, getNlocal()
1037	<
1036	>	do i = 1, nLocal
1037	>
1038		rfpot = 0.0_DP
1039		#ifdef IS_MPI
1040		me_i = atid_row(i)
1041		#else
1042		me_i = atid(i)
1043		#endif
1044	<	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
1045	<	if ( is_DP_i ) then
1046	<	call getElementProperty(atypes, me_i, "dipole_moment", mu_i)
1044	>
1045	>	if (PropertyMap(me_i)%is_DP) then
1046	>
1047	>	mu_i = PropertyMap(me_i)%dipole_moment
1048	>
1049		!! The reaction field needs to include a self contribution
1050		!! to the field:
1051	<	call accumulate_self_rf(i, mu_i, u_l)
1051	>	call accumulate_self_rf(i, mu_i, u_l)
1052		!! Get the reaction field contribution to the
1053		!! potential and torques:
1054		call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
#	Line 476 | Line 1062 | contains
1062		enddo
1063		endif
1064		endif
1065	<
1066	<
1065	>
1066	>
1067		#ifdef IS_MPI
1068	<
1068	>
1069		if (do_pot) then
1070		pot = pot + pot_local
1071		!! we assume the c code will do the allreduce to get the total potential
1072		!! we could do it right here if we needed to...
1073		endif
1074	<
1074	>
1075		if (do_stress) then
1076	<	call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1076	>	call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1077		mpi_comm_world,mpi_err)
1078		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1079		mpi_comm_world,mpi_err)
1080		endif
1081	<
1081	>
1082		#else
1083	<
1083	>
1084		if (do_stress) then
1085		tau = tau_Temp
1086		virial = virial_Temp
1087		endif
1088	<
1088	>
1089		#endif
1090
1091	<	write(*,*) 'T(1) = '
506	<	write(*,'(3es12.3)') t(1,1), t(1,2), t(1,3)
507	<	write(*,*)
508	<
509	<	write(*,*) 'T(2) = '
510	<	write(*,'(3es12.3)') t(2,1), t(2,2), t(2,3)
511	<	write(*,*)
512	<
1091	>
1092		end subroutine do_force_loop
1093
1094	<	subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
1094	>
1095	>	subroutine do_pair(i, j, rijsq, d, sw, do_pot, do_stress, &
1096	>	u_l, A, f, t, pot, vpair)
1097
1098	<	real( kind = dp ) :: pot
1099	<	real( kind = dp ), dimension(3,getNlocal()) :: u_l
1100	<	real (kind=dp), dimension(9,getNlocal()) :: A
1101	<	real (kind=dp), dimension(3,getNlocal()) :: f
1102	<	real (kind=dp), dimension(3,getNlocal()) :: t
1098	>	real( kind = dp ) :: pot, vpair, sw
1099	>	real( kind = dp ), dimension(nLocal)   :: mfact
1100	>	real( kind = dp ), dimension(3,nLocal) :: u_l
1101	>	real( kind = dp ), dimension(9,nLocal) :: A
1102	>	real( kind = dp ), dimension(3,nLocal) :: f
1103	>	real( kind = dp ), dimension(3,nLocal) :: t
1104
1105		logical, intent(inout) :: do_pot, do_stress
1106		integer, intent(in) :: i, j
1107	<	real ( kind = dp ), intent(inout)    :: rijsq
1107	>	real ( kind = dp ), intent(inout) :: rijsq
1108		real ( kind = dp )                :: r
1109		real ( kind = dp ), intent(inout) :: d(3)
528	–	logical :: is_LJ_i, is_LJ_j
529	–	logical :: is_DP_i, is_DP_j
530	–	logical :: is_GB_i, is_GB_j
531	–	logical :: is_Sticky_i, is_Sticky_j
1110		integer :: me_i, me_j
1111
1112		r = sqrt(rijsq)
1113
536	–	write(*,*) 'ul(1) = '
537	–	write(*,'(3es12.3)') u_l(1,1), u_l(1,2), u_l(1,3)
538	–	write(*,*)
539	–
540	–	write(*,*) 'ul(2) = '
541	–	write(*,'(3es12.3)') u_l(2,1), u_l(2,2), u_l(2,3)
542	–	write(*,*)
543	–
544	–
545	–	write(*,*) 'A(1) = '
546	–	write(*,'(3es12.3)') A(1,1), A(2,1), A(3,1)
547	–	write(*,'(3es12.3)') A(4,1), A(5,1), A(6,1)
548	–	write(*,'(3es12.3)') A(7,1), A(8,1), A(9,1)
549	–	write(*,*)
550	–	write(*,*) 'A(2) = '
551	–	write(*,'(3es12.3)') A(1,2), A(2,2), A(3,2)
552	–	write(*,'(3es12.3)') A(4,2), A(5,2), A(6,2)
553	–	write(*,'(3es12.3)') A(7,2), A(8,2), A(9,2)
554	–	write(*,*)
555	–
556	–
1114		#ifdef IS_MPI
1115		if (tagRow(i) .eq. tagColumn(j)) then
1116		write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
1117		endif
561	–
1118		me_i = atid_row(i)
1119		me_j = atid_col(j)
564	–
1120		#else
566	–
1121		me_i = atid(i)
1122		me_j = atid(j)
569	–
1123		#endif
1124	<
1125	<	if (FF_uses_LJ .and. SimUsesLJ()) then
1126	<	call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
1127	<	call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
1128	<
1129	<	if ( is_LJ_i .and. is_LJ_j ) &
1130	<	call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
1124	>
1125	>	if (FF_uses_LJ .and. SIM_uses_LJ) then
1126	>
1127	>	if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) then
1128	>	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, pot, f, do_pot, &
1129	>	do_stress)
1130	>	endif
1131	>
1132		endif
1133	<
1134	<	if (FF_uses_dipoles .and. SimUsesDipoles()) then
581	<	call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
582	<	call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
1133	>
1134	>	if (FF_uses_charges .and. SIM_uses_charges) then
1135
1136	<	if ( is_DP_i .and. is_DP_j ) then
1137	<
1138	<	call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
1136	>	if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then
1137	>	call do_charge_pair(i, j, d, r, rijsq, sw, vpair, pot, f, do_pot, &
1138	>	do_stress)
1139	>	endif
1140	>
1141	>	endif
1142	>
1143	>	if (FF_uses_dipoles .and. SIM_uses_dipoles) then
1144	>
1145	>	if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
1146	>	call do_dipole_pair(i, j, d, r, rijsq, sw, vpair, pot, u_l, f, t, &
1147		do_pot, do_stress)
1148	<	if (FF_uses_RF .and. SimUsesRF()) then
1148	>	if (FF_uses_RF .and. SIM_uses_RF) then
1149		call accumulate_rf(i, j, r, u_l)
1150		call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
1151	<	endif
592	<
1151	>	endif
1152		endif
1153	+
1154		endif
1155
1156	<	if (FF_uses_Sticky .and. SimUsesSticky()) then
1156	>	if (FF_uses_Sticky .and. SIM_uses_sticky) then
1157
1158	<	call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
1159	<	call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
600	<
601	<	if ( is_Sticky_i .and. is_Sticky_j ) then
602	<	call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
1158	>	if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
1159	>	call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, pot, A, f, t, &
1160		do_pot, do_stress)
1161		endif
1162	+
1163		endif
1164
1165
1166	<	if (FF_uses_GB .and. SimUsesGB()) then
609	<
610	<	call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
611	<	call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
1166	>	if (FF_uses_GB .and. SIM_uses_GB) then
1167
1168	<	if ( is_GB_i .and. is_GB_j ) then
1169	<	call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
1168	>	if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
1169	>	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, pot, u_l, f, t, &
1170		do_pot, do_stress)
1171		endif
617	–	endif
618	–
1172
1173	<
621	<	end subroutine do_pair
622	<
623	<
624	<	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
625	<
626	<	real (kind = dp), dimension(3) :: q_i
627	<	real (kind = dp), dimension(3) :: q_j
628	<	real ( kind = dp ), intent(out) :: r_sq
629	<	real( kind = dp ) :: d(3), scaled(3)
630	<	integer i
631	<
632	<	d(1:3) = q_j(1:3) - q_i(1:3)
633	<
634	<	! Wrap back into periodic box if necessary
635	<	if ( SimUsesPBC() ) then
1173	>	endif
1174
1175	<	if( .not.boxIsOrthorhombic ) then
1176	<	! calc the scaled coordinates.
1177	<
1178	<	scaled = matmul(HmatInv, d)
1179	<
642	<	! wrap the scaled coordinates
643	<
644	<	scaled = scaled  - anint(scaled)
645	<
646	<
647	<	! calc the wrapped real coordinates from the wrapped scaled
648	<	! coordinates
649	<
650	<	d = matmul(Hmat,scaled)
651	<
652	<	else
653	<	! calc the scaled coordinates.
654	<
655	<	do i = 1, 3
656	<	scaled(i) = d(i) * HmatInv(i,i)
657	<
658	<	! wrap the scaled coordinates
659	<
660	<	scaled(i) = scaled(i) - anint(scaled(i))
661	<
662	<	! calc the wrapped real coordinates from the wrapped scaled
663	<	! coordinates
664	<
665	<	d(i) = scaled(i)*Hmat(i,i)
666	<	enddo
1175	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1176	>
1177	>	if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
1178	>	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, pot, f, &
1179	>	do_pot, do_stress)
1180		endif
1181
1182		endif
1183
1184	<	r_sq = dot_product(d,d)
1185	<
1186	<	end subroutine get_interatomic_vector
1187	<
1188	<	subroutine check_initialization(error)
1189	<	integer, intent(out) :: error
1190	<
1191	<	error = 0
1192	<	! Make sure we are properly initialized.
1193	<	if (.not. do_forces_initialized) then
1194	<	error = -1
1195	<	return
1184	>	end subroutine do_pair
1185	>
1186	>	subroutine do_prepair(i, j, rijsq, d, rcijsq, dc, &
1187	>	do_pot, do_stress, u_l, A, f, t, pot)
1188	>	real( kind = dp ) :: pot
1189	>	real( kind = dp ), dimension(3,nLocal) :: u_l
1190	>	real (kind=dp), dimension(9,nLocal) :: A
1191	>	real (kind=dp), dimension(3,nLocal) :: f
1192	>	real (kind=dp), dimension(3,nLocal) :: t
1193	>
1194	>	logical, intent(inout) :: do_pot, do_stress
1195	>	integer, intent(in) :: i, j
1196	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1197	>	real ( kind = dp )                :: r, rc
1198	>	real ( kind = dp ), intent(inout) :: d(3), dc(3)
1199	>
1200	>	logical :: is_EAM_i, is_EAM_j
1201	>
1202	>	integer :: me_i, me_j
1203	>
1204	>
1205	>	r = sqrt(rijsq)
1206	>	if (SIM_uses_molecular_cutoffs) then
1207	>	rc = sqrt(rcijsq)
1208	>	else
1209	>	rc = r
1210		endif
1211	+
1212
1213		#ifdef IS_MPI
1214	<	if (.not. isMPISimSet()) then
1215	<	write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
1216	<	error = -1
1217	<	return
1218	<	endif
1214	>	if (tagRow(i) .eq. tagColumn(j)) then
1215	>	write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
1216	>	endif
1217	>
1218	>	me_i = atid_row(i)
1219	>	me_j = atid_col(j)
1220	>
1221	>	#else
1222	>
1223	>	me_i = atid(i)
1224	>	me_j = atid(j)
1225	>
1226		#endif
1227
1228	<	return
694	<	end subroutine check_initialization
1228	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1229
1230	<
1231	<	subroutine zero_work_arrays()
698	<
699	<	#ifdef IS_MPI
700	<
701	<	q_Row = 0.0_dp
702	<	q_Col = 0.0_dp
703	<
704	<	u_l_Row = 0.0_dp
705	<	u_l_Col = 0.0_dp
706	<
707	<	A_Row = 0.0_dp
708	<	A_Col = 0.0_dp
709	<
710	<	f_Row = 0.0_dp
711	<	f_Col = 0.0_dp
712	<	f_Temp = 0.0_dp
713	<
714	<	t_Row = 0.0_dp
715	<	t_Col = 0.0_dp
716	<	t_Temp = 0.0_dp
717	<
718	<	pot_Row = 0.0_dp
719	<	pot_Col = 0.0_dp
720	<	pot_Temp = 0.0_dp
1230	>	if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1231	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1232
1233	<	rf_Row = 0.0_dp
1234	<	rf_Col = 0.0_dp
1235	<	rf_Temp = 0.0_dp
1233	>	endif
1234	>
1235	>	end subroutine do_prepair
1236
726	–	#endif
1237
728	–	rf = 0.0_dp
729	–	tau_Temp = 0.0_dp
730	–	virial_Temp = 0.0_dp
731	–	end subroutine zero_work_arrays
732	–
733	–	function skipThisPair(atom1, atom2) result(skip_it)
734	–	integer, intent(in) :: atom1
735	–	integer, intent(in), optional :: atom2
736	–	logical :: skip_it
737	–	integer :: unique_id_1, unique_id_2
738	–	integer :: me_i,me_j
739	–	integer :: i
1238
1239	<	skip_it = .false.
1240	<
1241	<	!! there are a number of reasons to skip a pair or a particle
1242	<	!! mostly we do this to exclude atoms who are involved in short
1243	<	!! range interactions (bonds, bends, torsions), but we also need
1244	<	!! to exclude some overcounted interactions that result from
1245	<	!! the parallel decomposition
1246	<
1239	>
1240	>	subroutine do_preforce(nlocal,pot)
1241	>	integer :: nlocal
1242	>	real( kind = dp ) :: pot
1243	>
1244	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1245	>	call calc_EAM_preforce_Frho(nlocal,pot)
1246	>	endif
1247	>
1248	>
1249	>	end subroutine do_preforce
1250	>
1251	>
1252	>	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1253	>
1254	>	real (kind = dp), dimension(3) :: q_i
1255	>	real (kind = dp), dimension(3) :: q_j
1256	>	real ( kind = dp ), intent(out) :: r_sq
1257	>	real( kind = dp ) :: d(3), scaled(3)
1258	>	integer i
1259	>
1260	>	d(1:3) = q_j(1:3) - q_i(1:3)
1261	>
1262	>	! Wrap back into periodic box if necessary
1263	>	if ( SIM_uses_PBC ) then
1264	>
1265	>	if( .not.boxIsOrthorhombic ) then
1266	>	! calc the scaled coordinates.
1267	>
1268	>	scaled = matmul(HmatInv, d)
1269	>
1270	>	! wrap the scaled coordinates
1271	>
1272	>	scaled = scaled  - anint(scaled)
1273	>
1274	>
1275	>	! calc the wrapped real coordinates from the wrapped scaled
1276	>	! coordinates
1277	>
1278	>	d = matmul(Hmat,scaled)
1279	>
1280	>	else
1281	>	! calc the scaled coordinates.
1282	>
1283	>	do i = 1, 3
1284	>	scaled(i) = d(i) * HmatInv(i,i)
1285	>
1286	>	! wrap the scaled coordinates
1287	>
1288	>	scaled(i) = scaled(i) - anint(scaled(i))
1289	>
1290	>	! calc the wrapped real coordinates from the wrapped scaled
1291	>	! coordinates
1292	>
1293	>	d(i) = scaled(i)*Hmat(i,i)
1294	>	enddo
1295	>	endif
1296	>
1297	>	endif
1298	>
1299	>	r_sq = dot_product(d,d)
1300	>
1301	>	end subroutine get_interatomic_vector
1302	>
1303	>	subroutine zero_work_arrays()
1304	>
1305		#ifdef IS_MPI
1306	<	!! in MPI, we have to look up the unique IDs for each atom
1307	<	unique_id_1 = tagRow(atom1)
1306	>
1307	>	q_Row = 0.0_dp
1308	>	q_Col = 0.0_dp
1309	>
1310	>	q_group_Row = 0.0_dp
1311	>	q_group_Col = 0.0_dp
1312	>
1313	>	u_l_Row = 0.0_dp
1314	>	u_l_Col = 0.0_dp
1315	>
1316	>	A_Row = 0.0_dp
1317	>	A_Col = 0.0_dp
1318	>
1319	>	f_Row = 0.0_dp
1320	>	f_Col = 0.0_dp
1321	>	f_Temp = 0.0_dp
1322	>
1323	>	t_Row = 0.0_dp
1324	>	t_Col = 0.0_dp
1325	>	t_Temp = 0.0_dp
1326	>
1327	>	pot_Row = 0.0_dp
1328	>	pot_Col = 0.0_dp
1329	>	pot_Temp = 0.0_dp
1330	>
1331	>	rf_Row = 0.0_dp
1332	>	rf_Col = 0.0_dp
1333	>	rf_Temp = 0.0_dp
1334	>
1335	>	#endif
1336	>
1337	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1338	>	call clean_EAM()
1339	>	endif
1340	>
1341	>	rf = 0.0_dp
1342	>	tau_Temp = 0.0_dp
1343	>	virial_Temp = 0.0_dp
1344	>	end subroutine zero_work_arrays
1345	>
1346	>	function skipThisPair(atom1, atom2) result(skip_it)
1347	>	integer, intent(in) :: atom1
1348	>	integer, intent(in), optional :: atom2
1349	>	logical :: skip_it
1350	>	integer :: unique_id_1, unique_id_2
1351	>	integer :: me_i,me_j
1352	>	integer :: i
1353	>
1354	>	skip_it = .false.
1355	>
1356	>	!! there are a number of reasons to skip a pair or a particle
1357	>	!! mostly we do this to exclude atoms who are involved in short
1358	>	!! range interactions (bonds, bends, torsions), but we also need
1359	>	!! to exclude some overcounted interactions that result from
1360	>	!! the parallel decomposition
1361	>
1362	>	#ifdef IS_MPI
1363	>	!! in MPI, we have to look up the unique IDs for each atom
1364	>	unique_id_1 = tagRow(atom1)
1365		#else
1366	<	!! in the normal loop, the atom numbers are unique
1367	<	unique_id_1 = atom1
1366	>	!! in the normal loop, the atom numbers are unique
1367	>	unique_id_1 = atom1
1368		#endif
1369	<
1370	<	!! We were called with only one atom, so just check the global exclude
1371	<	!! list for this atom
1372	<	if (.not. present(atom2)) then
1373	<	do i = 1, nExcludes_global
1374	<	if (excludesGlobal(i) == unique_id_1) then
1375	<	skip_it = .true.
1376	<	return
1377	<	end if
1378	<	end do
1379	<	return
1380	<	end if
1381	<
1369	>
1370	>	!! We were called with only one atom, so just check the global exclude
1371	>	!! list for this atom
1372	>	if (.not. present(atom2)) then
1373	>	do i = 1, nExcludes_global
1374	>	if (excludesGlobal(i) == unique_id_1) then
1375	>	skip_it = .true.
1376	>	return
1377	>	end if
1378	>	end do
1379	>	return
1380	>	end if
1381	>
1382		#ifdef IS_MPI
1383	<	unique_id_2 = tagColumn(atom2)
1383	>	unique_id_2 = tagColumn(atom2)
1384		#else
1385	<	unique_id_2 = atom2
1385	>	unique_id_2 = atom2
1386		#endif
1387	<
1387	>
1388		#ifdef IS_MPI
1389	<	!! this situation should only arise in MPI simulations
1390	<	if (unique_id_1 == unique_id_2) then
1391	<	skip_it = .true.
1392	<	return
1393	<	end if
1394	<
1395	<	!! this prevents us from doing the pair on multiple processors
1396	<	if (unique_id_1 < unique_id_2) then
1397	<	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1398	<	skip_it = .true.
1399	<	return
1400	<	endif
1401	<	else
1402	<	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1403	<	skip_it = .true.
1404	<	return
1405	<	endif
1406	<	endif
1389	>	!! this situation should only arise in MPI simulations
1390	>	if (unique_id_1 == unique_id_2) then
1391	>	skip_it = .true.
1392	>	return
1393	>	end if
1394	>
1395	>	!! this prevents us from doing the pair on multiple processors
1396	>	if (unique_id_1 < unique_id_2) then
1397	>	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1398	>	skip_it = .true.
1399	>	return
1400	>	endif
1401	>	else
1402	>	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1403	>	skip_it = .true.
1404	>	return
1405	>	endif
1406	>	endif
1407		#endif
1408	+
1409	+	!! the rest of these situations can happen in all simulations:
1410	+	do i = 1, nExcludes_global
1411	+	if ((excludesGlobal(i) == unique_id_1) .or. &
1412	+	(excludesGlobal(i) == unique_id_2)) then
1413	+	skip_it = .true.
1414	+	return
1415	+	endif
1416	+	enddo
1417	+
1418	+	do i = 1, nExcludes_local
1419	+	if (excludesLocal(1,i) == unique_id_1) then
1420	+	if (excludesLocal(2,i) == unique_id_2) then
1421	+	skip_it = .true.
1422	+	return
1423	+	endif
1424	+	else
1425	+	if (excludesLocal(1,i) == unique_id_2) then
1426	+	if (excludesLocal(2,i) == unique_id_1) then
1427	+	skip_it = .true.
1428	+	return
1429	+	endif
1430	+	endif
1431	+	endif
1432	+	end do
1433	+
1434	+	return
1435	+	end function skipThisPair
1436
1437	<	!! the rest of these situations can happen in all simulations:
1438	<	do i = 1, nExcludes_global
1439	<	if ((excludesGlobal(i) == unique_id_1) .or. &
1440	<	(excludesGlobal(i) == unique_id_2)) then
1441	<	skip_it = .true.
1442	<	return
1443	<	endif
1444	<	enddo
1445	<
1446	<	do i = 1, nExcludes_local
1447	<	if (excludesLocal(1,i) == unique_id_1) then
1448	<	if (excludesLocal(2,i) == unique_id_2) then
1449	<	skip_it = .true.
1450	<	return
1451	<	endif
1452	<	else
1453	<	if (excludesLocal(1,i) == unique_id_2) then
1454	<	if (excludesLocal(2,i) == unique_id_1) then
1455	<	skip_it = .true.
1456	<	return
1457	<	endif
1458	<	endif
1459	<	endif
1460	<	end do
1461	<
821	<	return
822	<	end function skipThisPair
823	<
824	<	function FF_UsesDirectionalAtoms() result(doesit)
825	<	logical :: doesit
826	<	doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
827	<	FF_uses_GB .or. FF_uses_RF
828	<	end function FF_UsesDirectionalAtoms
829	<
830	<	function FF_RequiresPrepairCalc() result(doesit)
831	<	logical :: doesit
832	<	doesit = FF_uses_EAM
833	<	end function FF_RequiresPrepairCalc
834	<
835	<	function FF_RequiresPostpairCalc() result(doesit)
836	<	logical :: doesit
837	<	doesit = FF_uses_RF
838	<	end function FF_RequiresPostpairCalc
839	<
1437	>	function FF_UsesDirectionalAtoms() result(doesit)
1438	>	logical :: doesit
1439	>	doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
1440	>	FF_uses_GB .or. FF_uses_RF
1441	>	end function FF_UsesDirectionalAtoms
1442	>
1443	>	function FF_RequiresPrepairCalc() result(doesit)
1444	>	logical :: doesit
1445	>	doesit = FF_uses_EAM
1446	>	end function FF_RequiresPrepairCalc
1447	>
1448	>	function FF_RequiresPostpairCalc() result(doesit)
1449	>	logical :: doesit
1450	>	doesit = FF_uses_RF
1451	>	end function FF_RequiresPostpairCalc
1452	>
1453	>	#ifdef PROFILE
1454	>	function getforcetime() result(totalforcetime)
1455	>	real(kind=dp) :: totalforcetime
1456	>	totalforcetime = forcetime
1457	>	end function getforcetime
1458	>	#endif
1459	>
1460	>	!! This cleans componets of force arrays belonging only to fortran
1461	>
1462		end module do_Forces

Diff Legend

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