ViewVC Help

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

root/group/trunk/OOPSE-4/src/UseTheForce/doForces.F90

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 1628 by gezelter, Thu Oct 21 20:15:31 2004 UTC vs.
Revision 2259 by gezelter, Mon Jun 27 21:01:36 2005 UTC

#	Line 1 | Line 1
1	+	!!
2	+	!! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	+	!!
4	+	!! The University of Notre Dame grants you ("Licensee") a
5	+	!! non-exclusive, royalty free, license to use, modify and
6	+	!! redistribute this software in source and binary code form, provided
7	+	!! that the following conditions are met:
8	+	!!
9	+	!! 1. Acknowledgement of the program authors must be made in any
10	+	!!    publication of scientific results based in part on use of the
11	+	!!    program.  An acceptable form of acknowledgement is citation of
12	+	!!    the article in which the program was described (Matthew
13	+	!!    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	+	!!    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	+	!!    Parallel Simulation Engine for Molecular Dynamics,"
16	+	!!    J. Comput. Chem. 26, pp. 252-271 (2005))
17	+	!!
18	+	!! 2. Redistributions of source code must retain the above copyright
19	+	!!    notice, this list of conditions and the following disclaimer.
20	+	!!
21	+	!! 3. Redistributions in binary form must reproduce the above copyright
22	+	!!    notice, this list of conditions and the following disclaimer in the
23	+	!!    documentation and/or other materials provided with the
24	+	!!    distribution.
25	+	!!
26	+	!! This software is provided "AS IS," without a warranty of any
27	+	!! kind. All express or implied conditions, representations and
28	+	!! warranties, including any implied warranty of merchantability,
29	+	!! fitness for a particular purpose or non-infringement, are hereby
30	+	!! excluded.  The University of Notre Dame and its licensors shall not
31	+	!! be liable for any damages suffered by licensee as a result of
32	+	!! using, modifying or distributing the software or its
33	+	!! derivatives. In no event will the University of Notre Dame or its
34	+	!! licensors be liable for any lost revenue, profit or data, or for
35	+	!! direct, indirect, special, consequential, incidental or punitive
36	+	!! damages, however caused and regardless of the theory of liability,
37	+	!! arising out of the use of or inability to use software, even if the
38	+	!! University of Notre Dame has been advised of the possibility of
39	+	!! such damages.
40	+	!!
41	+
42		!! doForces.F90
43		!! module doForces
44		!! Calculates Long Range forces.
45
46		!! @author Charles F. Vardeman II
47		!! @author Matthew Meineke
48	<	!! @version $Id: doForces.F90,v 1.2 2004-10-21 20:15:22 gezelter Exp $, $Date: 2004-10-21 20:15:22 $, $Name: not supported by cvs2svn $, $Revision: 1.2 $
48	>	!! @version $Id: doForces.F90,v 1.20 2005-06-27 21:01:30 gezelter Exp $, $Date: 2005-06-27 21:01:30 $, $Name: not supported by cvs2svn $, $Revision: 1.20 $
49
50	+
51		module doForces
52		use force_globals
53		use simulation
#	Line 14 | Line 56 | module doForces
56		use switcheroo
57		use neighborLists
58		use lj
59	<	use sticky_pair
60	<	use dipole_dipole
19	<	use charge_charge
59	>	use sticky
60	>	use electrostatic_module
61		use reaction_field
62		use gb_pair
63	+	use shapes
64		use vector_class
65		use eam
66		use status
#	Line 31 | Line 73 | module doForces
73
74		#define __FORTRAN90
75		#include "UseTheForce/fSwitchingFunction.h"
76	+	#include "UseTheForce/DarkSide/fInteractionMap.h"
77
78		INTEGER, PARAMETER:: PREPAIR_LOOP = 1
79		INTEGER, PARAMETER:: PAIR_LOOP    = 2
#	Line 40 | Line 83 | module doForces
83		logical, save :: haveSIMvariables = .false.
84		logical, save :: havePropertyMap = .false.
85		logical, save :: haveSaneForceField = .false.
86	<	logical, save :: FF_uses_LJ
87	<	logical, save :: FF_uses_sticky
88	<	logical, save :: FF_uses_charges
89	<	logical, save :: FF_uses_dipoles
90	<	logical, save :: FF_uses_RF
91	<	logical, save :: FF_uses_GB
86	>
87	>	logical, save :: FF_uses_DirectionalAtoms
88	>	logical, save :: FF_uses_LennardJones
89	>	logical, save :: FF_uses_Electrostatics
90	>	logical, save :: FF_uses_Charges
91	>	logical, save :: FF_uses_Dipoles
92	>	logical, save :: FF_uses_Quadrupoles
93	>	logical, save :: FF_uses_Sticky
94	>	logical, save :: FF_uses_StickyPower
95	>	logical, save :: FF_uses_GayBerne
96		logical, save :: FF_uses_EAM
97	<	logical, save :: SIM_uses_LJ
98	<	logical, save :: SIM_uses_sticky
99	<	logical, save :: SIM_uses_charges
100	<	logical, save :: SIM_uses_dipoles
101	<	logical, save :: SIM_uses_RF
102	<	logical, save :: SIM_uses_GB
97	>	logical, save :: FF_uses_Shapes
98	>	logical, save :: FF_uses_FLARB
99	>	logical, save :: FF_uses_RF
100	>
101	>	logical, save :: SIM_uses_DirectionalAtoms
102	>	logical, save :: SIM_uses_LennardJones
103	>	logical, save :: SIM_uses_Electrostatics
104	>	logical, save :: SIM_uses_Charges
105	>	logical, save :: SIM_uses_Dipoles
106	>	logical, save :: SIM_uses_Quadrupoles
107	>	logical, save :: SIM_uses_Sticky
108	>	logical, save :: SIM_uses_StickyPower
109	>	logical, save :: SIM_uses_GayBerne
110		logical, save :: SIM_uses_EAM
111	+	logical, save :: SIM_uses_Shapes
112	+	logical, save :: SIM_uses_FLARB
113	+	logical, save :: SIM_uses_RF
114		logical, save :: SIM_requires_postpair_calc
115		logical, save :: SIM_requires_prepair_calc
59	–	logical, save :: SIM_uses_directional_atoms
116		logical, save :: SIM_uses_PBC
117		logical, save :: SIM_uses_molecular_cutoffs
118
#	Line 74 | Line 130 | module doForces
130		#endif
131
132		type :: Properties
133	<	logical :: is_lj     = .false.
134	<	logical :: is_sticky = .false.
135	<	logical :: is_dp     = .false.
136	<	logical :: is_gb     = .false.
137	<	logical :: is_eam    = .false.
138	<	logical :: is_charge = .false.
139	<	real(kind=DP) :: charge = 0.0_DP
140	<	real(kind=DP) :: dipole_moment = 0.0_DP
133	>	logical :: is_Directional   = .false.
134	>	logical :: is_LennardJones  = .false.
135	>	logical :: is_Electrostatic = .false.
136	>	logical :: is_Charge        = .false.
137	>	logical :: is_Dipole        = .false.
138	>	logical :: is_Quadrupole    = .false.
139	>	logical :: is_Sticky        = .false.
140	>	logical :: is_StickyPower   = .false.
141	>	logical :: is_GayBerne      = .false.
142	>	logical :: is_EAM           = .false.
143	>	logical :: is_Shape         = .false.
144	>	logical :: is_FLARB         = .false.
145		end type Properties
146
147		type(Properties), dimension(:),allocatable :: PropertyMap
#	Line 89 | Line 149 | contains
149		contains
150
151		subroutine setRlistDF( this_rlist )
152	<
152	>
153		real(kind=dp) :: this_rlist
154
155		rlist = this_rlist
156		rlistsq = rlist * rlist
157	<
157	>
158		haveRlist = .true.
159
160	<	end subroutine setRlistDF
160	>	end subroutine setRlistDF
161
162		subroutine createPropertyMap(status)
163		integer :: nAtypes
#	Line 114 | Line 174 | contains
174		status = -1
175		return
176		end if
177	<
177	>
178		if (.not. allocated(PropertyMap)) then
179		allocate(PropertyMap(nAtypes))
180		endif
181
182		do i = 1, nAtypes
183	<	call getElementProperty(atypes, i, "is_LJ", thisProperty)
184	<	PropertyMap(i)%is_LJ = thisProperty
183	>	call getElementProperty(atypes, i, "is_Directional", thisProperty)
184	>	PropertyMap(i)%is_Directional = thisProperty
185
186	+	call getElementProperty(atypes, i, "is_LennardJones", thisProperty)
187	+	PropertyMap(i)%is_LennardJones = thisProperty
188	+
189	+	call getElementProperty(atypes, i, "is_Electrostatic", thisProperty)
190	+	PropertyMap(i)%is_Electrostatic = thisProperty
191	+
192		call getElementProperty(atypes, i, "is_Charge", thisProperty)
193		PropertyMap(i)%is_Charge = thisProperty
128	–
129	–	if (thisProperty) then
130	–	call getElementProperty(atypes, i, "charge", thisDPproperty)
131	–	PropertyMap(i)%charge = thisDPproperty
132	–	endif
194
195	<	call getElementProperty(atypes, i, "is_DP", thisProperty)
196	<	PropertyMap(i)%is_DP = thisProperty
195	>	call getElementProperty(atypes, i, "is_Dipole", thisProperty)
196	>	PropertyMap(i)%is_Dipole = thisProperty
197
198	<	if (thisProperty) then
199	<	call getElementProperty(atypes, i, "dipole_moment", thisDPproperty)
139	<	PropertyMap(i)%dipole_moment = thisDPproperty
140	<	endif
198	>	call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
199	>	PropertyMap(i)%is_Quadrupole = thisProperty
200
201		call getElementProperty(atypes, i, "is_Sticky", thisProperty)
202		PropertyMap(i)%is_Sticky = thisProperty
203	<	call getElementProperty(atypes, i, "is_GB", thisProperty)
204	<	PropertyMap(i)%is_GB = thisProperty
203	>
204	>	call getElementProperty(atypes, i, "is_StickyPower", thisProperty)
205	>	PropertyMap(i)%is_StickyPower = thisProperty
206	>
207	>	call getElementProperty(atypes, i, "is_GayBerne", thisProperty)
208	>	PropertyMap(i)%is_GayBerne = thisProperty
209	>
210		call getElementProperty(atypes, i, "is_EAM", thisProperty)
211		PropertyMap(i)%is_EAM = thisProperty
212	+
213	+	call getElementProperty(atypes, i, "is_Shape", thisProperty)
214	+	PropertyMap(i)%is_Shape = thisProperty
215	+
216	+	call getElementProperty(atypes, i, "is_FLARB", thisProperty)
217	+	PropertyMap(i)%is_FLARB = thisProperty
218		end do
219
220		havePropertyMap = .true.
#	Line 152 | Line 222 | contains
222		end subroutine createPropertyMap
223
224		subroutine setSimVariables()
225	<	SIM_uses_LJ = SimUsesLJ()
226	<	SIM_uses_sticky = SimUsesSticky()
227	<	SIM_uses_charges = SimUsesCharges()
228	<	SIM_uses_dipoles = SimUsesDipoles()
229	<	SIM_uses_RF = SimUsesRF()
230	<	SIM_uses_GB = SimUsesGB()
225	>	SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
226	>	SIM_uses_LennardJones = SimUsesLennardJones()
227	>	SIM_uses_Electrostatics = SimUsesElectrostatics()
228	>	SIM_uses_Charges = SimUsesCharges()
229	>	SIM_uses_Dipoles = SimUsesDipoles()
230	>	SIM_uses_Sticky = SimUsesSticky()
231	>	SIM_uses_StickyPower = SimUsesStickyPower()
232	>	SIM_uses_GayBerne = SimUsesGayBerne()
233		SIM_uses_EAM = SimUsesEAM()
234	+	SIM_uses_Shapes = SimUsesShapes()
235	+	SIM_uses_FLARB = SimUsesFLARB()
236	+	SIM_uses_RF = SimUsesRF()
237		SIM_requires_postpair_calc = SimRequiresPostpairCalc()
238		SIM_requires_prepair_calc = SimRequiresPrepairCalc()
164	–	SIM_uses_directional_atoms = SimUsesDirectionalAtoms()
239		SIM_uses_PBC = SimUsesPBC()
166	–	!SIM_uses_molecular_cutoffs = SimUsesMolecularCutoffs()
240
241		haveSIMvariables = .true.
242
#	Line 176 | Line 249 | contains
249		integer :: myStatus
250
251		error = 0
252	<
252	>
253		if (.not. havePropertyMap) then
254
255		myStatus = 0
#	Line 221 | Line 294 | contains
294		#endif
295		return
296		end subroutine doReadyCheck
224	–
297
298	+
299		subroutine init_FF(use_RF_c, thisStat)
300
301		logical, intent(in) :: use_RF_c
#	Line 237 | Line 310 | contains
310
311		!! Fortran's version of a cast:
312		FF_uses_RF = use_RF_c
313	<
313	>
314		!! init_FF is called *after* all of the atom types have been
315		!! defined in atype_module using the new_atype subroutine.
316		!!
317		!! this will scan through the known atypes and figure out what
318		!! interactions are used by the force field.
319	<
320	<	FF_uses_LJ = .false.
321	<	FF_uses_sticky = .false.
322	<	FF_uses_charges = .false.
323	<	FF_uses_dipoles = .false.
324	<	FF_uses_GB = .false.
319	>
320	>	FF_uses_DirectionalAtoms = .false.
321	>	FF_uses_LennardJones = .false.
322	>	FF_uses_Electrostatics = .false.
323	>	FF_uses_Charges = .false.
324	>	FF_uses_Dipoles = .false.
325	>	FF_uses_Sticky = .false.
326	>	FF_uses_StickyPower = .false.
327	>	FF_uses_GayBerne = .false.
328		FF_uses_EAM = .false.
329	<
330	<	call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList)
331	<	if (nMatches .gt. 0) FF_uses_LJ = .true.
332	<
333	<	call getMatchingElementList(atypes, "is_Charge", .true., nMatches, MatchList)
334	<	if (nMatches .gt. 0) FF_uses_charges = .true.
335	<
336	<	call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
337	<	if (nMatches .gt. 0) FF_uses_dipoles = .true.
338	<
329	>	FF_uses_Shapes = .false.
330	>	FF_uses_FLARB = .false.
331	>
332	>	call getMatchingElementList(atypes, "is_Directional", .true., &
333	>	nMatches, MatchList)
334	>	if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
335	>
336	>	call getMatchingElementList(atypes, "is_LennardJones", .true., &
337	>	nMatches, MatchList)
338	>	if (nMatches .gt. 0) FF_uses_LennardJones = .true.
339	>
340	>	call getMatchingElementList(atypes, "is_Electrostatic", .true., &
341	>	nMatches, MatchList)
342	>	if (nMatches .gt. 0) then
343	>	FF_uses_Electrostatics = .true.
344	>	endif
345	>
346	>	call getMatchingElementList(atypes, "is_Charge", .true., &
347	>	nMatches, MatchList)
348	>	if (nMatches .gt. 0) then
349	>	FF_uses_Charges = .true.
350	>	FF_uses_Electrostatics = .true.
351	>	endif
352	>
353	>	call getMatchingElementList(atypes, "is_Dipole", .true., &
354	>	nMatches, MatchList)
355	>	if (nMatches .gt. 0) then
356	>	FF_uses_Dipoles = .true.
357	>	FF_uses_Electrostatics = .true.
358	>	FF_uses_DirectionalAtoms = .true.
359	>	endif
360	>
361	>	call getMatchingElementList(atypes, "is_Quadrupole", .true., &
362	>	nMatches, MatchList)
363	>	if (nMatches .gt. 0) then
364	>	FF_uses_Quadrupoles = .true.
365	>	FF_uses_Electrostatics = .true.
366	>	FF_uses_DirectionalAtoms = .true.
367	>	endif
368	>
369		call getMatchingElementList(atypes, "is_Sticky", .true., nMatches, &
370		MatchList)
371	<	if (nMatches .gt. 0) FF_uses_Sticky = .true.
371	>	if (nMatches .gt. 0) then
372	>	FF_uses_Sticky = .true.
373	>	FF_uses_DirectionalAtoms = .true.
374	>	endif
375	>
376	>	call getMatchingElementList(atypes, "is_StickyPower", .true., nMatches, &
377	>	MatchList)
378	>	if (nMatches .gt. 0) then
379	>	FF_uses_StickyPower = .true.
380	>	FF_uses_DirectionalAtoms = .true.
381	>	endif
382
383	<	call getMatchingElementList(atypes, "is_GB", .true., nMatches, MatchList)
384	<	if (nMatches .gt. 0) FF_uses_GB = .true.
385	<
383	>	call getMatchingElementList(atypes, "is_GayBerne", .true., &
384	>	nMatches, MatchList)
385	>	if (nMatches .gt. 0) then
386	>	FF_uses_GayBerne = .true.
387	>	FF_uses_DirectionalAtoms = .true.
388	>	endif
389	>
390		call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
391		if (nMatches .gt. 0) FF_uses_EAM = .true.
392	<
392	>
393	>	call getMatchingElementList(atypes, "is_Shape", .true., &
394	>	nMatches, MatchList)
395	>	if (nMatches .gt. 0) then
396	>	FF_uses_Shapes = .true.
397	>	FF_uses_DirectionalAtoms = .true.
398	>	endif
399	>
400	>	call getMatchingElementList(atypes, "is_FLARB", .true., &
401	>	nMatches, MatchList)
402	>	if (nMatches .gt. 0) FF_uses_FLARB = .true.
403	>
404		!! Assume sanity (for the sake of argument)
405		haveSaneForceField = .true.
406	<
406	>
407		!! check to make sure the FF_uses_RF setting makes sense
408	<
408	>
409		if (FF_uses_dipoles) then
410		if (FF_uses_RF) then
411		dielect = getDielect()
#	Line 287 | Line 418 | contains
418		haveSaneForceField = .false.
419		return
420		endif
290	–	endif
291	–
292	–	if (FF_uses_sticky) then
293	–	call check_sticky_FF(my_status)
294	–	if (my_status /= 0) then
295	–	thisStat = -1
296	–	haveSaneForceField = .false.
297	–	return
298	–	end if
421		endif
422
423	+	!sticky module does not contain check_sticky_FF anymore
424	+	!if (FF_uses_sticky) then
425	+	!   call check_sticky_FF(my_status)
426	+	!   if (my_status /= 0) then
427	+	!      thisStat = -1
428	+	!      haveSaneForceField = .false.
429	+	!      return
430	+	!   end if
431	+	!endif
432	+
433		if (FF_uses_EAM) then
434	<	call init_EAM_FF(my_status)
434	>	call init_EAM_FF(my_status)
435		if (my_status /= 0) then
436		write(default_error, *) "init_EAM_FF returned a bad status"
437		thisStat = -1
#	Line 308 | Line 440 | contains
440		end if
441		endif
442
443	<	if (FF_uses_GB) then
443	>	if (FF_uses_GayBerne) then
444		call check_gb_pair_FF(my_status)
445		if (my_status .ne. 0) then
446		thisStat = -1
#	Line 317 | Line 449 | contains
449		endif
450		endif
451
452	<	if (FF_uses_GB .and. FF_uses_LJ) then
452	>	if (FF_uses_GayBerne .and. FF_uses_LennardJones) then
453		endif
454
455		if (.not. haveNeighborList) then
#	Line 329 | Line 461 | contains
461		return
462		endif
463		haveNeighborList = .true.
464	<	endif
465	<
464	>	endif
465	>
466		end subroutine init_FF
335	–
467
468	+
469		!! Does force loop over i,j pairs. Calls do_pair to calculates forces.
470		!------------------------------------------------------------->
471	<	subroutine do_force_loop(q, q_group, A, u_l, f, t, tau, pot, &
471	>	subroutine do_force_loop(q, q_group, A, eFrame, f, t, tau, pot, &
472		do_pot_c, do_stress_c, error)
473		!! Position array provided by C, dimensioned by getNlocal
474		real ( kind = dp ), dimension(3, nLocal) :: q
#	Line 345 | Line 477 | contains
477		!! Rotation Matrix for each long range particle in simulation.
478		real( kind = dp), dimension(9, nLocal) :: A
479		!! Unit vectors for dipoles (lab frame)
480	<	real( kind = dp ), dimension(3,nLocal) :: u_l
480	>	real( kind = dp ), dimension(9,nLocal) :: eFrame
481		!! Force array provided by C, dimensioned by getNlocal
482		real ( kind = dp ), dimension(3,nLocal) :: f
483		!! Torsion array provided by C, dimensioned by getNlocal
#	Line 385 | Line 517 | contains
517		integer :: loopStart, loopEnd, loop
518
519		real(kind=dp) :: listSkin = 1.0
520	<
520	>
521		!! initialize local variables
522	<
522	>
523		#ifdef IS_MPI
524		pot_local = 0.0_dp
525		nAtomsInRow   = getNatomsInRow(plan_atom_row)
#	Line 397 | Line 529 | contains
529		#else
530		natoms = nlocal
531		#endif
532	<
532	>
533		call doReadyCheck(localError)
534		if ( localError .ne. 0 ) then
535		call handleError("do_force_loop", "Not Initialized")
#	Line 405 | Line 537 | contains
537		return
538		end if
539		call zero_work_arrays()
540	<
540	>
541		do_pot = do_pot_c
542		do_stress = do_stress_c
543	<
543	>
544		! Gather all information needed by all force loops:
545	<
545	>
546		#ifdef IS_MPI
547	<
547	>
548		call gather(q, q_Row, plan_atom_row_3d)
549		call gather(q, q_Col, plan_atom_col_3d)
550
551		call gather(q_group, q_group_Row, plan_group_row_3d)
552		call gather(q_group, q_group_Col, plan_group_col_3d)
553	<
554	<	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
555	<	call gather(u_l, u_l_Row, plan_atom_row_3d)
556	<	call gather(u_l, u_l_Col, plan_atom_col_3d)
557	<
553	>
554	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
555	>	call gather(eFrame, eFrame_Row, plan_atom_row_rotation)
556	>	call gather(eFrame, eFrame_Col, plan_atom_col_rotation)
557	>
558		call gather(A, A_Row, plan_atom_row_rotation)
559		call gather(A, A_Col, plan_atom_col_rotation)
560		endif
561	<
561	>
562		#endif
563	<
563	>
564		!! Begin force loop timing:
565		#ifdef PROFILE
566		call cpu_time(forceTimeInitial)
567		nloops = nloops + 1
568		#endif
569	<
569	>
570		loopEnd = PAIR_LOOP
571		if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
572		loopStart = PREPAIR_LOOP
#	Line 449 | Line 581 | contains
581		if (loop .eq. loopStart) then
582		#ifdef IS_MPI
583		call checkNeighborList(nGroupsInRow, q_group_row, listSkin, &
584	<	update_nlist)
584	>	update_nlist)
585		#else
586		call checkNeighborList(nGroups, q_group, listSkin, &
587	<	update_nlist)
587	>	update_nlist)
588		#endif
589		endif
590	<
590	>
591		if (update_nlist) then
592		!! save current configuration and construct neighbor list
593		#ifdef IS_MPI
#	Line 466 | Line 598 | contains
598		neighborListSize = size(list)
599		nlist = 0
600		endif
601	<
601	>
602		istart = 1
603		#ifdef IS_MPI
604		iend = nGroupsInRow
#	Line 476 | Line 608 | contains
608		outer: do i = istart, iend
609
610		if (update_nlist) point(i) = nlist + 1
611	<
611	>
612		n_in_i = groupStartRow(i+1) - groupStartRow(i)
613	<
613	>
614		if (update_nlist) then
615		#ifdef IS_MPI
616		jstart = 1
#	Line 493 | Line 625 | contains
625		! make sure group i has neighbors
626		if (jstart .gt. jend) cycle outer
627		endif
628	<
628	>
629		do jnab = jstart, jend
630		if (update_nlist) then
631		j = jnab
#	Line 512 | Line 644 | contains
644		if (rgrpsq < rlistsq) then
645		if (update_nlist) then
646		nlist = nlist + 1
647	<
647	>
648		if (nlist > neighborListSize) then
649		#ifdef IS_MPI
650		call expandNeighborList(nGroupsInRow, listerror)
#	Line 526 | Line 658 | contains
658		end if
659		neighborListSize = size(list)
660		endif
661	<
661	>
662		list(nlist) = j
663		endif
664	<
664	>
665		if (loop .eq. PAIR_LOOP) then
666		vij = 0.0d0
667		fij(1:3) = 0.0d0
668		endif
669	<
669	>
670		call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
671		in_switching_region)
672	<
672	>
673		n_in_j = groupStartCol(j+1) - groupStartCol(j)
674	<
674	>
675		do ia = groupStartRow(i), groupStartRow(i+1)-1
676	<
676	>
677		atom1 = groupListRow(ia)
678	<
678	>
679		inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
680	<
680	>
681		atom2 = groupListCol(jb)
682	<
682	>
683		if (skipThisPair(atom1, atom2)) cycle inner
684
685		if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
#	Line 567 | Line 699 | contains
699		#ifdef IS_MPI
700		call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
701		rgrpsq, d_grp, do_pot, do_stress, &
702	<	u_l, A, f, t, pot_local)
702	>	eFrame, A, f, t, pot_local)
703		#else
704		call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
705		rgrpsq, d_grp, do_pot, do_stress, &
706	<	u_l, A, f, t, pot)
706	>	eFrame, A, f, t, pot)
707		#endif
708		else
709		#ifdef IS_MPI
710		call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
711		do_pot, &
712	<	u_l, A, f, t, pot_local, vpair, fpair)
712	>	eFrame, A, f, t, pot_local, vpair, fpair)
713		#else
714		call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
715		do_pot,  &
716	<	u_l, A, f, t, pot, vpair, fpair)
716	>	eFrame, A, f, t, pot, vpair, fpair)
717		#endif
718
719		vij = vij + vpair
#	Line 589 | Line 721 | contains
721		endif
722		enddo inner
723		enddo
724	<
724	>
725		if (loop .eq. PAIR_LOOP) then
726		if (in_switching_region) then
727		swderiv = vij*dswdr/rgrp
728		fij(1) = fij(1) + swderiv*d_grp(1)
729		fij(2) = fij(2) + swderiv*d_grp(2)
730		fij(3) = fij(3) + swderiv*d_grp(3)
731	<
731	>
732		do ia=groupStartRow(i), groupStartRow(i+1)-1
733		atom1=groupListRow(ia)
734		mf = mfactRow(atom1)
#	Line 610 | Line 742 | contains
742		f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
743		#endif
744		enddo
745	<
745	>
746		do jb=groupStartCol(j), groupStartCol(j+1)-1
747		atom2=groupListCol(jb)
748		mf = mfactCol(atom2)
#	Line 625 | Line 757 | contains
757		#endif
758		enddo
759		endif
760	<
760	>
761		if (do_stress) call add_stress_tensor(d_grp, fij)
762		endif
763		end if
764		enddo
765		enddo outer
766	<
766	>
767		if (update_nlist) then
768		#ifdef IS_MPI
769		point(nGroupsInRow + 1) = nlist + 1
#	Line 645 | Line 777 | contains
777		update_nlist = .false.
778		endif
779		endif
780	<
780	>
781		if (loop .eq. PREPAIR_LOOP) then
782		call do_preforce(nlocal, pot)
783		endif
784	<
784	>
785		enddo
786	<
786	>
787		!! Do timing
788		#ifdef PROFILE
789		call cpu_time(forceTimeFinal)
790		forceTime = forceTime + forceTimeFinal - forceTimeInitial
791		#endif
792	<
792	>
793		#ifdef IS_MPI
794		!!distribute forces
795	<
795	>
796		f_temp = 0.0_dp
797		call scatter(f_Row,f_temp,plan_atom_row_3d)
798		do i = 1,nlocal
799		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
800		end do
801	<
801	>
802		f_temp = 0.0_dp
803		call scatter(f_Col,f_temp,plan_atom_col_3d)
804		do i = 1,nlocal
805		f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
806		end do
807	<
808	<	if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
807	>
808	>	if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
809		t_temp = 0.0_dp
810		call scatter(t_Row,t_temp,plan_atom_row_3d)
811		do i = 1,nlocal
#	Line 681 | Line 813 | contains
813		end do
814		t_temp = 0.0_dp
815		call scatter(t_Col,t_temp,plan_atom_col_3d)
816	<
816	>
817		do i = 1,nlocal
818		t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
819		end do
820		endif
821	<
821	>
822		if (do_pot) then
823		! scatter/gather pot_row into the members of my column
824		call scatter(pot_Row, pot_Temp, plan_atom_row)
825	<
825	>
826		! scatter/gather pot_local into all other procs
827		! add resultant to get total pot
828		do i = 1, nlocal
829		pot_local = pot_local + pot_Temp(i)
830		enddo
831	<
831	>
832		pot_Temp = 0.0_DP
833	<
833	>
834		call scatter(pot_Col, pot_Temp, plan_atom_col)
835		do i = 1, nlocal
836		pot_local = pot_local + pot_Temp(i)
837		enddo
838	<
838	>
839		endif
840		#endif
841	<
841	>
842		if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
843	<
843	>
844		if (FF_uses_RF .and. SIM_uses_RF) then
845	<
845	>
846		#ifdef IS_MPI
847		call scatter(rf_Row,rf,plan_atom_row_3d)
848		call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
#	Line 718 | Line 850 | contains
850		rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
851		end do
852		#endif
853	<
853	>
854		do i = 1, nLocal
855	<
855	>
856		rfpot = 0.0_DP
857		#ifdef IS_MPI
858		me_i = atid_row(i)
859		#else
860		me_i = atid(i)
861		#endif
862	<
863	<	if (PropertyMap(me_i)%is_DP) then
864	<
865	<	mu_i = PropertyMap(me_i)%dipole_moment
866	<
862	>
863	>	if (PropertyMap(me_i)%is_Dipole) then
864	>
865	>	mu_i = getDipoleMoment(me_i)
866	>
867		!! The reaction field needs to include a self contribution
868		!! to the field:
869	<	call accumulate_self_rf(i, mu_i, u_l)
869	>	call accumulate_self_rf(i, mu_i, eFrame)
870		!! Get the reaction field contribution to the
871		!! potential and torques:
872	<	call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
872	>	call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
873		#ifdef IS_MPI
874		pot_local = pot_local + rfpot
875		#else
876		pot = pot + rfpot
877	<
877	>
878		#endif
879	<	endif
879	>	endif
880		enddo
881		endif
882		endif
883	<
884	<
883	>
884	>
885		#ifdef IS_MPI
886	<
886	>
887		if (do_pot) then
888		pot = pot + pot_local
889		!! we assume the c code will do the allreduce to get the total potential
890		!! we could do it right here if we needed to...
891		endif
892	<
892	>
893		if (do_stress) then
894		call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
895		mpi_comm_world,mpi_err)
896		call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
897		mpi_comm_world,mpi_err)
898		endif
899	<
899	>
900		#else
901	<
901	>
902		if (do_stress) then
903		tau = tau_Temp
904		virial = virial_Temp
905		endif
906	<
906	>
907		#endif
908	<
908	>
909		end subroutine do_force_loop
910	<
910	>
911		subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
912	<	u_l, A, f, t, pot, vpair, fpair)
912	>	eFrame, A, f, t, pot, vpair, fpair)
913
914		real( kind = dp ) :: pot, vpair, sw
915		real( kind = dp ), dimension(3) :: fpair
916		real( kind = dp ), dimension(nLocal)   :: mfact
917	<	real( kind = dp ), dimension(3,nLocal) :: u_l
917	>	real( kind = dp ), dimension(9,nLocal) :: eFrame
918		real( kind = dp ), dimension(9,nLocal) :: A
919		real( kind = dp ), dimension(3,nLocal) :: f
920		real( kind = dp ), dimension(3,nLocal) :: t
#	Line 792 | Line 924 | contains
924		real ( kind = dp ), intent(inout) :: rijsq
925		real ( kind = dp )                :: r
926		real ( kind = dp ), intent(inout) :: d(3)
927	+	real ( kind = dp ) :: ebalance
928		integer :: me_i, me_j
929
930	+	integer :: iMap
931	+
932		r = sqrt(rijsq)
933		vpair = 0.0d0
934		fpair(1:3) = 0.0d0
#	Line 805 | Line 940 | contains
940		me_i = atid(i)
941		me_j = atid(j)
942		#endif
808	–
809	–	if (FF_uses_LJ .and. SIM_uses_LJ) then
810	–
811	–	if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) then
812	–	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
813	–	endif
814	–
815	–	endif
816	–
817	–	if (FF_uses_charges .and. SIM_uses_charges) then
818	–
819	–	if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then
820	–	call do_charge_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
821	–	endif
822	–
823	–	endif
824	–
825	–	if (FF_uses_dipoles .and. SIM_uses_dipoles) then
826	–
827	–	if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
828	–	call do_dipole_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, u_l, f, t, &
829	–	do_pot)
830	–	if (FF_uses_RF .and. SIM_uses_RF) then
831	–	call accumulate_rf(i, j, r, u_l, sw)
832	–	call rf_correct_forces(i, j, d, r, u_l, sw, f, fpair)
833	–	endif
834	–	endif
943
944	+	iMap = InteractionMap(me_i, me_j)%InteractionHash
945	+
946	+	if ( iand(iMap, LJ_PAIR).ne.0 ) then
947	+	call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
948		endif
949
950	<	if (FF_uses_Sticky .and. SIM_uses_sticky) then
950	>	if ( iand(iMap, ELECTROSTATIC_PAIR).ne.0 ) then
951	>	call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
952	>	pot, eFrame, f, t, do_pot)
953
954	<	if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
955	<	call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, A, f, t, &
956	<	do_pot)
954	>	if (FF_uses_dipoles .and. SIM_uses_dipoles) then
955	>	if ( PropertyMap(me_i)%is_Dipole .and. &
956	>	PropertyMap(me_j)%is_Dipole) then
957	>	if (FF_uses_RF .and. SIM_uses_RF) then
958	>	call accumulate_rf(i, j, r, eFrame, sw)
959	>	call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
960	>	endif
961	>	endif
962		endif
963	+	endif
964
965	+	if ( iand(iMap, STICKY_PAIR).ne.0 ) then
966	+	call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
967	+	pot, A, f, t, do_pot)
968		endif
969
970	+	if ( iand(iMap, STICKYPOWER_PAIR).ne.0 ) then
971	+	call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
972	+	pot, A, f, t, do_pot)
973	+	endif
974
975	<	if (FF_uses_GB .and. SIM_uses_GB) then
976	<
977	<	if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
978	<	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, u_l, f, t, &
979	<	do_pot)
980	<	endif
975	>	if ( iand(iMap, GAYBERNE_PAIR).ne.0 ) then
976	>	call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
977	>	pot, A, f, t, do_pot)
978	>	endif
979	>
980	>	if ( iand(iMap, GAYBERNE_LJ).ne.0 ) then
981	>	call do_gblj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
982	>	pot, A, f, t, do_pot)
983	>	endif
984
985	+	if ( iand(iMap, EAM_PAIR).ne.0 ) then
986	+	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
987	+	do_pot)
988		endif
989	<
990	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
991	<
992	<	if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
860	<	call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
861	<	do_pot)
862	<	endif
863	<
989	>
990	>	if ( iand(iMap, SHAPE_PAIR).ne.0 ) then
991	>	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
992	>	pot, A, f, t, do_pot)
993		endif
994	+
995	+	if ( iand(iMap, SHAPE_LJ).ne.0 ) then
996	+	call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
997	+	pot, A, f, t, do_pot)
998	+	endif
999
1000		end subroutine do_pair
1001
1002		subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1003	<	do_pot, do_stress, u_l, A, f, t, pot)
1003	>	do_pot, do_stress, eFrame, A, f, t, pot)
1004
1005	<	real( kind = dp ) :: pot, sw
1006	<	real( kind = dp ), dimension(3,nLocal) :: u_l
1007	<	real (kind=dp), dimension(9,nLocal) :: A
1008	<	real (kind=dp), dimension(3,nLocal) :: f
1009	<	real (kind=dp), dimension(3,nLocal) :: t
876	<
877	<	logical, intent(inout) :: do_pot, do_stress
878	<	integer, intent(in) :: i, j
879	<	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
880	<	real ( kind = dp )                :: r, rc
881	<	real ( kind = dp ), intent(inout) :: d(3), dc(3)
882	<
883	<	logical :: is_EAM_i, is_EAM_j
884	<
885	<	integer :: me_i, me_j
886	<
1005	>	real( kind = dp ) :: pot, sw
1006	>	real( kind = dp ), dimension(9,nLocal) :: eFrame
1007	>	real (kind=dp), dimension(9,nLocal) :: A
1008	>	real (kind=dp), dimension(3,nLocal) :: f
1009	>	real (kind=dp), dimension(3,nLocal) :: t
1010
1011	<	r = sqrt(rijsq)
1012	<	if (SIM_uses_molecular_cutoffs) then
1013	<	rc = sqrt(rcijsq)
1014	<	else
1015	<	rc = r
893	<	endif
894	<
1011	>	logical, intent(inout) :: do_pot, do_stress
1012	>	integer, intent(in) :: i, j
1013	>	real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
1014	>	real ( kind = dp )                :: r, rc
1015	>	real ( kind = dp ), intent(inout) :: d(3), dc(3)
1016
1017	+	integer :: me_i, me_j, iMap
1018	+
1019		#ifdef IS_MPI
1020	<	me_i = atid_row(i)
1021	<	me_j = atid_col(j)
1020	>	me_i = atid_row(i)
1021	>	me_j = atid_col(j)
1022		#else
1023	<	me_i = atid(i)
1024	<	me_j = atid(j)
1023	>	me_i = atid(i)
1024	>	me_j = atid(j)
1025		#endif
1026	<
1027	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1028	<
1029	<	if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1030	<	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1031	<
1032	<	endif
1033	<
1034	<	end subroutine do_prepair
1035	<
1036	<
1037	<	subroutine do_preforce(nlocal,pot)
1038	<	integer :: nlocal
1039	<	real( kind = dp ) :: pot
1040	<
1041	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1042	<	call calc_EAM_preforce_Frho(nlocal,pot)
1043	<	endif
1044	<
1045	<
1046	<	end subroutine do_preforce
1047	<
1048	<
1049	<	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1050	<
1051	<	real (kind = dp), dimension(3) :: q_i
1052	<	real (kind = dp), dimension(3) :: q_j
1053	<	real ( kind = dp ), intent(out) :: r_sq
1054	<	real( kind = dp ) :: d(3), scaled(3)
1055	<	integer i
1056	<
1057	<	d(1:3) = q_j(1:3) - q_i(1:3)
1058	<
1059	<	! Wrap back into periodic box if necessary
1060	<	if ( SIM_uses_PBC ) then
1061	<
1062	<	if( .not.boxIsOrthorhombic ) then
1063	<	! calc the scaled coordinates.
1064	<
1065	<	scaled = matmul(HmatInv, d)
1066	<
1067	<	! wrap the scaled coordinates
1068	<
1069	<	scaled = scaled  - anint(scaled)
1070	<
1071	<
1072	<	! calc the wrapped real coordinates from the wrapped scaled
1073	<	! coordinates
1074	<
1075	<	d = matmul(Hmat,scaled)
1076	<
1077	<	else
1078	<	! calc the scaled coordinates.
1079	<
1080	<	do i = 1, 3
1081	<	scaled(i) = d(i) * HmatInv(i,i)
1082	<
1083	<	! wrap the scaled coordinates
1084	<
1085	<	scaled(i) = scaled(i) - anint(scaled(i))
1086	<
1087	<	! calc the wrapped real coordinates from the wrapped scaled
1088	<	! coordinates
1089	<
1090	<	d(i) = scaled(i)*Hmat(i,i)
1091	<	enddo
1092	<	endif
1093	<
1094	<	endif
1095	<
1096	<	r_sq = dot_product(d,d)
1097	<
1098	<	end subroutine get_interatomic_vector
1099	<
1100	<	subroutine zero_work_arrays()
978	<
1026	>
1027	>	iMap = InteractionMap(me_i, me_j)%InteractionHash
1028	>
1029	>	if ( iand(iMap, EAM_PAIR).ne.0 ) then
1030	>	call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1031	>	endif
1032	>
1033	>	end subroutine do_prepair
1034	>
1035	>
1036	>	subroutine do_preforce(nlocal,pot)
1037	>	integer :: nlocal
1038	>	real( kind = dp ) :: pot
1039	>
1040	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1041	>	call calc_EAM_preforce_Frho(nlocal,pot)
1042	>	endif
1043	>
1044	>
1045	>	end subroutine do_preforce
1046	>
1047	>
1048	>	subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1049	>
1050	>	real (kind = dp), dimension(3) :: q_i
1051	>	real (kind = dp), dimension(3) :: q_j
1052	>	real ( kind = dp ), intent(out) :: r_sq
1053	>	real( kind = dp ) :: d(3), scaled(3)
1054	>	integer i
1055	>
1056	>	d(1:3) = q_j(1:3) - q_i(1:3)
1057	>
1058	>	! Wrap back into periodic box if necessary
1059	>	if ( SIM_uses_PBC ) then
1060	>
1061	>	if( .not.boxIsOrthorhombic ) then
1062	>	! calc the scaled coordinates.
1063	>
1064	>	scaled = matmul(HmatInv, d)
1065	>
1066	>	! wrap the scaled coordinates
1067	>
1068	>	scaled = scaled  - anint(scaled)
1069	>
1070	>
1071	>	! calc the wrapped real coordinates from the wrapped scaled
1072	>	! coordinates
1073	>
1074	>	d = matmul(Hmat,scaled)
1075	>
1076	>	else
1077	>	! calc the scaled coordinates.
1078	>
1079	>	do i = 1, 3
1080	>	scaled(i) = d(i) * HmatInv(i,i)
1081	>
1082	>	! wrap the scaled coordinates
1083	>
1084	>	scaled(i) = scaled(i) - anint(scaled(i))
1085	>
1086	>	! calc the wrapped real coordinates from the wrapped scaled
1087	>	! coordinates
1088	>
1089	>	d(i) = scaled(i)*Hmat(i,i)
1090	>	enddo
1091	>	endif
1092	>
1093	>	endif
1094	>
1095	>	r_sq = dot_product(d,d)
1096	>
1097	>	end subroutine get_interatomic_vector
1098	>
1099	>	subroutine zero_work_arrays()
1100	>
1101		#ifdef IS_MPI
980	–
981	–	q_Row = 0.0_dp
982	–	q_Col = 0.0_dp
1102
1103	<	q_group_Row = 0.0_dp
1104	<	q_group_Col = 0.0_dp
1105	<
1106	<	u_l_Row = 0.0_dp
1107	<	u_l_Col = 0.0_dp
1108	<
1109	<	A_Row = 0.0_dp
1110	<	A_Col = 0.0_dp
1111	<
1112	<	f_Row = 0.0_dp
1113	<	f_Col = 0.0_dp
1114	<	f_Temp = 0.0_dp
1115	<
1116	<	t_Row = 0.0_dp
1117	<	t_Col = 0.0_dp
1118	<	t_Temp = 0.0_dp
1119	<
1120	<	pot_Row = 0.0_dp
1121	<	pot_Col = 0.0_dp
1122	<	pot_Temp = 0.0_dp
1123	<
1124	<	rf_Row = 0.0_dp
1125	<	rf_Col = 0.0_dp
1126	<	rf_Temp = 0.0_dp
1127	<
1103	>	q_Row = 0.0_dp
1104	>	q_Col = 0.0_dp
1105	>
1106	>	q_group_Row = 0.0_dp
1107	>	q_group_Col = 0.0_dp
1108	>
1109	>	eFrame_Row = 0.0_dp
1110	>	eFrame_Col = 0.0_dp
1111	>
1112	>	A_Row = 0.0_dp
1113	>	A_Col = 0.0_dp
1114	>
1115	>	f_Row = 0.0_dp
1116	>	f_Col = 0.0_dp
1117	>	f_Temp = 0.0_dp
1118	>
1119	>	t_Row = 0.0_dp
1120	>	t_Col = 0.0_dp
1121	>	t_Temp = 0.0_dp
1122	>
1123	>	pot_Row = 0.0_dp
1124	>	pot_Col = 0.0_dp
1125	>	pot_Temp = 0.0_dp
1126	>
1127	>	rf_Row = 0.0_dp
1128	>	rf_Col = 0.0_dp
1129	>	rf_Temp = 0.0_dp
1130	>
1131		#endif
1132	<
1133	<	if (FF_uses_EAM .and. SIM_uses_EAM) then
1134	<	call clean_EAM()
1135	<	endif
1136	<
1137	<	rf = 0.0_dp
1138	<	tau_Temp = 0.0_dp
1139	<	virial_Temp = 0.0_dp
1140	<	end subroutine zero_work_arrays
1141	<
1142	<	function skipThisPair(atom1, atom2) result(skip_it)
1143	<	integer, intent(in) :: atom1
1144	<	integer, intent(in), optional :: atom2
1145	<	logical :: skip_it
1146	<	integer :: unique_id_1, unique_id_2
1147	<	integer :: me_i,me_j
1148	<	integer :: i
1149	<
1150	<	skip_it = .false.
1151	<
1152	<	!! there are a number of reasons to skip a pair or a particle
1153	<	!! mostly we do this to exclude atoms who are involved in short
1154	<	!! range interactions (bonds, bends, torsions), but we also need
1155	<	!! to exclude some overcounted interactions that result from
1156	<	!! the parallel decomposition
1157	<
1158	<	#ifdef IS_MPI
1159	<	!! in MPI, we have to look up the unique IDs for each atom
1160	<	unique_id_1 = AtomRowToGlobal(atom1)
1161	<	#else
1162	<	!! in the normal loop, the atom numbers are unique
1163	<	unique_id_1 = atom1
1132	>
1133	>	if (FF_uses_EAM .and. SIM_uses_EAM) then
1134	>	call clean_EAM()
1135	>	endif
1136	>
1137	>	rf = 0.0_dp
1138	>	tau_Temp = 0.0_dp
1139	>	virial_Temp = 0.0_dp
1140	>	end subroutine zero_work_arrays
1141	>
1142	>	function skipThisPair(atom1, atom2) result(skip_it)
1143	>	integer, intent(in) :: atom1
1144	>	integer, intent(in), optional :: atom2
1145	>	logical :: skip_it
1146	>	integer :: unique_id_1, unique_id_2
1147	>	integer :: me_i,me_j
1148	>	integer :: i
1149	>
1150	>	skip_it = .false.
1151	>
1152	>	!! there are a number of reasons to skip a pair or a particle
1153	>	!! mostly we do this to exclude atoms who are involved in short
1154	>	!! range interactions (bonds, bends, torsions), but we also need
1155	>	!! to exclude some overcounted interactions that result from
1156	>	!! the parallel decomposition
1157	>
1158	>	#ifdef IS_MPI
1159	>	!! in MPI, we have to look up the unique IDs for each atom
1160	>	unique_id_1 = AtomRowToGlobal(atom1)
1161	>	#else
1162	>	!! in the normal loop, the atom numbers are unique
1163	>	unique_id_1 = atom1
1164		#endif
1165	<
1166	<	!! We were called with only one atom, so just check the global exclude
1167	<	!! list for this atom
1168	<	if (.not. present(atom2)) then
1169	<	do i = 1, nExcludes_global
1170	<	if (excludesGlobal(i) == unique_id_1) then
1171	<	skip_it = .true.
1172	<	return
1173	<	end if
1174	<	end do
1175	<	return
1176	<	end if
1177	<
1165	>
1166	>	!! We were called with only one atom, so just check the global exclude
1167	>	!! list for this atom
1168	>	if (.not. present(atom2)) then
1169	>	do i = 1, nExcludes_global
1170	>	if (excludesGlobal(i) == unique_id_1) then
1171	>	skip_it = .true.
1172	>	return
1173	>	end if
1174	>	end do
1175	>	return
1176	>	end if
1177	>
1178		#ifdef IS_MPI
1179	<	unique_id_2 = AtomColToGlobal(atom2)
1179	>	unique_id_2 = AtomColToGlobal(atom2)
1180		#else
1181	<	unique_id_2 = atom2
1181	>	unique_id_2 = atom2
1182		#endif
1183	<
1183	>
1184		#ifdef IS_MPI
1185	<	!! this situation should only arise in MPI simulations
1186	<	if (unique_id_1 == unique_id_2) then
1187	<	skip_it = .true.
1188	<	return
1189	<	end if
1190	<
1191	<	!! this prevents us from doing the pair on multiple processors
1192	<	if (unique_id_1 < unique_id_2) then
1193	<	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1194	<	skip_it = .true.
1195	<	return
1196	<	endif
1197	<	else
1198	<	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1199	<	skip_it = .true.
1200	<	return
1201	<	endif
1202	<	endif
1185	>	!! this situation should only arise in MPI simulations
1186	>	if (unique_id_1 == unique_id_2) then
1187	>	skip_it = .true.
1188	>	return
1189	>	end if
1190	>
1191	>	!! this prevents us from doing the pair on multiple processors
1192	>	if (unique_id_1 < unique_id_2) then
1193	>	if (mod(unique_id_1 + unique_id_2,2) == 0) then
1194	>	skip_it = .true.
1195	>	return
1196	>	endif
1197	>	else
1198	>	if (mod(unique_id_1 + unique_id_2,2) == 1) then
1199	>	skip_it = .true.
1200	>	return
1201	>	endif
1202	>	endif
1203		#endif
1204	<
1205	<	!! the rest of these situations can happen in all simulations:
1206	<	do i = 1, nExcludes_global
1207	<	if ((excludesGlobal(i) == unique_id_1) .or. &
1208	<	(excludesGlobal(i) == unique_id_2)) then
1209	<	skip_it = .true.
1210	<	return
1211	<	endif
1212	<	enddo
1213	<
1214	<	do i = 1, nSkipsForAtom(atom1)
1215	<	if (skipsForAtom(atom1, i) .eq. unique_id_2) then
1216	<	skip_it = .true.
1217	<	return
1218	<	endif
1219	<	end do
1220	<
1221	<	return
1222	<	end function skipThisPair
1223	<
1224	<	function FF_UsesDirectionalAtoms() result(doesit)
1225	<	logical :: doesit
1226	<	doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
1227	<	FF_uses_GB .or. FF_uses_RF
1228	<	end function FF_UsesDirectionalAtoms
1229	<
1230	<	function FF_RequiresPrepairCalc() result(doesit)
1231	<	logical :: doesit
1232	<	doesit = FF_uses_EAM
1233	<	end function FF_RequiresPrepairCalc
1234	<
1235	<	function FF_RequiresPostpairCalc() result(doesit)
1236	<	logical :: doesit
1237	<	doesit = FF_uses_RF
1238	<	end function FF_RequiresPostpairCalc
1239	<
1204	>
1205	>	!! the rest of these situations can happen in all simulations:
1206	>	do i = 1, nExcludes_global
1207	>	if ((excludesGlobal(i) == unique_id_1) .or. &
1208	>	(excludesGlobal(i) == unique_id_2)) then
1209	>	skip_it = .true.
1210	>	return
1211	>	endif
1212	>	enddo
1213	>
1214	>	do i = 1, nSkipsForAtom(atom1)
1215	>	if (skipsForAtom(atom1, i) .eq. unique_id_2) then
1216	>	skip_it = .true.
1217	>	return
1218	>	endif
1219	>	end do
1220	>
1221	>	return
1222	>	end function skipThisPair
1223	>
1224	>	function FF_UsesDirectionalAtoms() result(doesit)
1225	>	logical :: doesit
1226	>	doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1227	>	FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1228	>	FF_uses_StickyPower .or. FF_uses_GayBerne .or. FF_uses_Shapes
1229	>	end function FF_UsesDirectionalAtoms
1230	>
1231	>	function FF_RequiresPrepairCalc() result(doesit)
1232	>	logical :: doesit
1233	>	doesit = FF_uses_EAM
1234	>	end function FF_RequiresPrepairCalc
1235	>
1236	>	function FF_RequiresPostpairCalc() result(doesit)
1237	>	logical :: doesit
1238	>	doesit = FF_uses_RF
1239	>	end function FF_RequiresPostpairCalc
1240	>
1241		#ifdef PROFILE
1242	<	function getforcetime() result(totalforcetime)
1243	<	real(kind=dp) :: totalforcetime
1244	<	totalforcetime = forcetime
1245	<	end function getforcetime
1242	>	function getforcetime() result(totalforcetime)
1243	>	real(kind=dp) :: totalforcetime
1244	>	totalforcetime = forcetime
1245	>	end function getforcetime
1246		#endif
1124	–
1125	–	!! This cleans componets of force arrays belonging only to fortran
1247
1248	<	subroutine add_stress_tensor(dpair, fpair)
1128	<
1129	<	real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1130	<
1131	<	! because the d vector is the rj - ri vector, and
1132	<	! because fx, fy, fz are the force on atom i, we need a
1133	<	! negative sign here:
1134	<
1135	<	tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1136	<	tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1137	<	tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1138	<	tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1139	<	tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1140	<	tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1141	<	tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1142	<	tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1143	<	tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1144	<
1145	<	virial_Temp = virial_Temp + &
1146	<	(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1147	<
1148	<	end subroutine add_stress_tensor
1149	<
1150	<	end module doForces
1248	>	!! This cleans componets of force arrays belonging only to fortran
1249
1250	<	!! Interfaces for C programs to module....
1250	>	subroutine add_stress_tensor(dpair, fpair)
1251
1252	<	subroutine initFortranFF(use_RF_c, thisStat)
1155	<	use doForces, ONLY: init_FF
1156	<	logical, intent(in) :: use_RF_c
1252	>	real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1253
1254	<	integer, intent(out) :: thisStat
1255	<	call init_FF(use_RF_c, thisStat)
1254	>	! because the d vector is the rj - ri vector, and
1255	>	! because fx, fy, fz are the force on atom i, we need a
1256	>	! negative sign here:
1257
1258	<	end subroutine initFortranFF
1258	>	tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1259	>	tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1260	>	tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1261	>	tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1262	>	tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1263	>	tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1264	>	tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1265	>	tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1266	>	tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1267
1268	<	subroutine doForceloop(q, q_group, A, u_l, f, t, tau, pot, &
1269	<	do_pot_c, do_stress_c, error)
1165	<
1166	<	use definitions, ONLY: dp
1167	<	use simulation
1168	<	use doForces, ONLY: do_force_loop
1169	<	!! Position array provided by C, dimensioned by getNlocal
1170	<	real ( kind = dp ), dimension(3, nLocal) :: q
1171	<	!! molecular center-of-mass position array
1172	<	real ( kind = dp ), dimension(3, nGroups) :: q_group
1173	<	!! Rotation Matrix for each long range particle in simulation.
1174	<	real( kind = dp), dimension(9, nLocal) :: A
1175	<	!! Unit vectors for dipoles (lab frame)
1176	<	real( kind = dp ), dimension(3,nLocal) :: u_l
1177	<	!! Force array provided by C, dimensioned by getNlocal
1178	<	real ( kind = dp ), dimension(3,nLocal) :: f
1179	<	!! Torsion array provided by C, dimensioned by getNlocal
1180	<	real( kind = dp ), dimension(3,nLocal) :: t
1268	>	virial_Temp = virial_Temp + &
1269	>	(tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1270
1271	<	!! Stress Tensor
1272	<	real( kind = dp), dimension(9) :: tau
1273	<	real ( kind = dp ) :: pot
1185	<	logical ( kind = 2) :: do_pot_c, do_stress_c
1186	<	integer :: error
1187	<
1188	<	call do_force_loop(q, q_group, A, u_l, f, t, tau, pot, &
1189	<	do_pot_c, do_stress_c, error)
1190	<
1191	<	end subroutine doForceloop
1271	>	end subroutine add_stress_tensor
1272	>
1273	>	end module doForces

Diff Legend

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