ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/doForces.F90
Revision: 937
Committed: Sun Apr 16 02:51:16 2006 UTC (19 years ago) by chrisfen
File size: 49997 byte(s)
Log Message:
added a cubic spline to switcheroo

File Contents

# User Rev Content
1 gezelter 246 !!
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 gezelter 117 !! doForces.F90
43     !! module doForces
44     !! Calculates Long Range forces.
45    
46     !! @author Charles F. Vardeman II
47     !! @author Matthew Meineke
48 chrisfen 937 !! @version $Id: doForces.F90,v 1.77 2006-04-16 02:51:16 chrisfen Exp $, $Date: 2006-04-16 02:51:16 $, $Name: not supported by cvs2svn $, $Revision: 1.77 $
49 gezelter 117
50 gezelter 246
51 gezelter 117 module doForces
52     use force_globals
53     use simulation
54     use definitions
55     use atype_module
56     use switcheroo
57     use neighborLists
58     use lj
59 gezelter 246 use sticky
60 gezelter 401 use electrostatic_module
61 gezelter 676 use gayberne
62 chrisfen 143 use shapes
63 gezelter 117 use vector_class
64     use eam
65 chuckv 733 use suttonchen
66 gezelter 117 use status
67 chrisfen 937 use interpolation
68 gezelter 117 #ifdef IS_MPI
69     use mpiSimulation
70     #endif
71    
72     implicit none
73     PRIVATE
74    
75     #define __FORTRAN90
76     #include "UseTheForce/fSwitchingFunction.h"
77 gezelter 574 #include "UseTheForce/fCutoffPolicy.h"
78 gezelter 560 #include "UseTheForce/DarkSide/fInteractionMap.h"
79 chrisfen 611 #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
80 gezelter 117
81 gezelter 574
82 gezelter 117 INTEGER, PARAMETER:: PREPAIR_LOOP = 1
83     INTEGER, PARAMETER:: PAIR_LOOP = 2
84    
85     logical, save :: haveNeighborList = .false.
86     logical, save :: haveSIMvariables = .false.
87     logical, save :: haveSaneForceField = .false.
88 gezelter 571 logical, save :: haveInteractionHash = .false.
89     logical, save :: haveGtypeCutoffMap = .false.
90 chrisfen 618 logical, save :: haveDefaultCutoffs = .false.
91 gezelter 762 logical, save :: haveSkinThickness = .false.
92     logical, save :: haveElectrostaticSummationMethod = .false.
93     logical, save :: haveCutoffPolicy = .false.
94     logical, save :: VisitCutoffsAfterComputing = .false.
95 chrisfen 532
96 gezelter 141 logical, save :: FF_uses_DirectionalAtoms
97 gezelter 401 logical, save :: FF_uses_Dipoles
98 gezelter 141 logical, save :: FF_uses_GayBerne
99     logical, save :: FF_uses_EAM
100 chuckv 733 logical, save :: FF_uses_SC
101     logical, save :: FF_uses_MEAM
102    
103 gezelter 141
104     logical, save :: SIM_uses_DirectionalAtoms
105     logical, save :: SIM_uses_EAM
106 chuckv 733 logical, save :: SIM_uses_SC
107     logical, save :: SIM_uses_MEAM
108 gezelter 117 logical, save :: SIM_requires_postpair_calc
109     logical, save :: SIM_requires_prepair_calc
110     logical, save :: SIM_uses_PBC
111    
112 chrisfen 607 integer, save :: electrostaticSummationMethod
113 gezelter 762 integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
114 chrisfen 580
115 gezelter 762 real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut
116     real(kind=dp), save :: skinThickness
117     logical, save :: defaultDoShift
118    
119 gezelter 117 public :: init_FF
120 gezelter 762 public :: setCutoffs
121     public :: cWasLame
122     public :: setElectrostaticMethod
123     public :: setCutoffPolicy
124     public :: setSkinThickness
125 gezelter 117 public :: do_force_loop
126    
127     #ifdef PROFILE
128     public :: getforcetime
129     real, save :: forceTime = 0
130     real :: forceTimeInitial, forceTimeFinal
131     integer :: nLoops
132     #endif
133 chuckv 561
134 gezelter 571 !! Variables for cutoff mapping and interaction mapping
135     ! Bit hash to determine pair-pair interactions.
136     integer, dimension(:,:), allocatable :: InteractionHash
137     real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
138 chuckv 651 real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
139     real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
140    
141     integer, dimension(:), allocatable, target :: groupToGtypeRow
142     integer, dimension(:), pointer :: groupToGtypeCol => null()
143    
144     real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
145     real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
146 gezelter 571 type ::gtypeCutoffs
147     real(kind=dp) :: rcut
148     real(kind=dp) :: rcutsq
149     real(kind=dp) :: rlistsq
150     end type gtypeCutoffs
151     type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
152 gezelter 574
153 gezelter 117 contains
154    
155 gezelter 762 subroutine createInteractionHash()
156 chuckv 561 integer :: nAtypes
157     integer :: i
158     integer :: j
159 gezelter 571 integer :: iHash
160 tim 568 !! Test Types
161 chuckv 561 logical :: i_is_LJ
162     logical :: i_is_Elect
163     logical :: i_is_Sticky
164     logical :: i_is_StickyP
165     logical :: i_is_GB
166     logical :: i_is_EAM
167     logical :: i_is_Shape
168 chuckv 733 logical :: i_is_SC
169     logical :: i_is_MEAM
170 chuckv 561 logical :: j_is_LJ
171     logical :: j_is_Elect
172     logical :: j_is_Sticky
173     logical :: j_is_StickyP
174     logical :: j_is_GB
175     logical :: j_is_EAM
176     logical :: j_is_Shape
177 chuckv 733 logical :: j_is_SC
178     logical :: j_is_MEAM
179 gezelter 576 real(kind=dp) :: myRcut
180    
181 chuckv 561 if (.not. associated(atypes)) then
182 gezelter 762 call handleError("doForces", "atypes was not present before call of createInteractionHash!")
183 chuckv 561 return
184     endif
185    
186     nAtypes = getSize(atypes)
187    
188     if (nAtypes == 0) then
189 gezelter 762 call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
190 chuckv 561 return
191     end if
192 chrisfen 532
193 chuckv 570 if (.not. allocated(InteractionHash)) then
194     allocate(InteractionHash(nAtypes,nAtypes))
195 chuckv 655 else
196     deallocate(InteractionHash)
197     allocate(InteractionHash(nAtypes,nAtypes))
198 chuckv 561 endif
199 gezelter 571
200     if (.not. allocated(atypeMaxCutoff)) then
201     allocate(atypeMaxCutoff(nAtypes))
202 chuckv 655 else
203     deallocate(atypeMaxCutoff)
204     allocate(atypeMaxCutoff(nAtypes))
205 gezelter 571 endif
206 chuckv 561
207     do i = 1, nAtypes
208     call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
209     call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
210     call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
211     call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
212     call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
213     call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
214     call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
215 chuckv 733 call getElementProperty(atypes, i, "is_SC", i_is_SC)
216     call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM)
217 gezelter 117
218 chuckv 561 do j = i, nAtypes
219 chrisfen 532
220 chuckv 561 iHash = 0
221     myRcut = 0.0_dp
222 gezelter 117
223 chuckv 561 call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
224     call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect)
225     call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky)
226     call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP)
227     call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
228     call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
229     call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
230 chuckv 733 call getElementProperty(atypes, j, "is_SC", j_is_SC)
231     call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM)
232 gezelter 117
233 chuckv 561 if (i_is_LJ .and. j_is_LJ) then
234 gezelter 562 iHash = ior(iHash, LJ_PAIR)
235     endif
236    
237     if (i_is_Elect .and. j_is_Elect) then
238     iHash = ior(iHash, ELECTROSTATIC_PAIR)
239     endif
240    
241     if (i_is_Sticky .and. j_is_Sticky) then
242     iHash = ior(iHash, STICKY_PAIR)
243     endif
244 chuckv 561
245 gezelter 562 if (i_is_StickyP .and. j_is_StickyP) then
246     iHash = ior(iHash, STICKYPOWER_PAIR)
247     endif
248 chuckv 561
249 gezelter 562 if (i_is_EAM .and. j_is_EAM) then
250     iHash = ior(iHash, EAM_PAIR)
251 chuckv 561 endif
252    
253 chuckv 733 if (i_is_SC .and. j_is_SC) then
254     iHash = ior(iHash, SC_PAIR)
255     endif
256    
257 chuckv 561 if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
258     if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
259     if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
260    
261     if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
262     if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
263     if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
264    
265    
266 chuckv 570 InteractionHash(i,j) = iHash
267     InteractionHash(j,i) = iHash
268 chuckv 561
269     end do
270    
271     end do
272 tim 568
273 gezelter 571 haveInteractionHash = .true.
274     end subroutine createInteractionHash
275 chuckv 561
276 gezelter 762 subroutine createGtypeCutoffMap()
277 gezelter 569
278 gezelter 574 logical :: i_is_LJ
279     logical :: i_is_Elect
280     logical :: i_is_Sticky
281     logical :: i_is_StickyP
282     logical :: i_is_GB
283     logical :: i_is_EAM
284     logical :: i_is_Shape
285 chuckv 831 logical :: i_is_SC
286 gezelter 587 logical :: GtypeFound
287 chuckv 561
288 gezelter 576 integer :: myStatus, nAtypes, i, j, istart, iend, jstart, jend
289 chuckv 652 integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
290 chuckv 589 integer :: nGroupsInRow
291 chuckv 651 integer :: nGroupsInCol
292     integer :: nGroupTypesRow,nGroupTypesCol
293 gezelter 762 real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
294 gezelter 576 real(kind=dp) :: biggestAtypeCutoff
295 gezelter 571
296     if (.not. haveInteractionHash) then
297 gezelter 762 call createInteractionHash()
298 chuckv 567 endif
299 chuckv 589 #ifdef IS_MPI
300     nGroupsInRow = getNgroupsInRow(plan_group_row)
301 chuckv 651 nGroupsInCol = getNgroupsInCol(plan_group_col)
302 chuckv 589 #endif
303 chuckv 563 nAtypes = getSize(atypes)
304 chuckv 599 ! Set all of the initial cutoffs to zero.
305     atypeMaxCutoff = 0.0_dp
306 gezelter 571 do i = 1, nAtypes
307 gezelter 582 if (SimHasAtype(i)) then
308 gezelter 575 call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
309     call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
310     call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
311     call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
312     call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
313     call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
314     call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
315 chuckv 831 call getElementProperty(atypes, i, "is_SC", i_is_SC)
316 chuckv 599
317 chrisfen 618 if (haveDefaultCutoffs) then
318     atypeMaxCutoff(i) = defaultRcut
319     else
320     if (i_is_LJ) then
321     thisRcut = getSigma(i) * 2.5_dp
322     if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
323     endif
324     if (i_is_Elect) then
325     thisRcut = defaultRcut
326     if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
327     endif
328     if (i_is_Sticky) then
329     thisRcut = getStickyCut(i)
330     if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
331     endif
332     if (i_is_StickyP) then
333     thisRcut = getStickyPowerCut(i)
334     if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
335     endif
336     if (i_is_GB) then
337     thisRcut = getGayBerneCut(i)
338     if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
339     endif
340     if (i_is_EAM) then
341     thisRcut = getEAMCut(i)
342     if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
343     endif
344     if (i_is_Shape) then
345     thisRcut = getShapeCut(i)
346     if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
347     endif
348 chuckv 831 if (i_is_SC) then
349     thisRcut = getSCCut(i)
350     if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
351     endif
352 gezelter 575 endif
353 gezelter 762
354 gezelter 575 if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
355     biggestAtypeCutoff = atypeMaxCutoff(i)
356     endif
357 chrisfen 618
358 gezelter 574 endif
359 gezelter 575 enddo
360 gezelter 581
361 gezelter 575 istart = 1
362 chuckv 651 jstart = 1
363 gezelter 575 #ifdef IS_MPI
364     iend = nGroupsInRow
365 chuckv 651 jend = nGroupsInCol
366 gezelter 575 #else
367     iend = nGroups
368 chuckv 651 jend = nGroups
369 gezelter 575 #endif
370 gezelter 582
371 gezelter 581 !! allocate the groupToGtype and gtypeMaxCutoff here.
372 chuckv 651 if(.not.allocated(groupToGtypeRow)) then
373     ! allocate(groupToGtype(iend))
374     allocate(groupToGtypeRow(iend))
375     else
376     deallocate(groupToGtypeRow)
377     allocate(groupToGtypeRow(iend))
378 chuckv 583 endif
379 chuckv 651 if(.not.allocated(groupMaxCutoffRow)) then
380     allocate(groupMaxCutoffRow(iend))
381     else
382     deallocate(groupMaxCutoffRow)
383     allocate(groupMaxCutoffRow(iend))
384     end if
385    
386     if(.not.allocated(gtypeMaxCutoffRow)) then
387     allocate(gtypeMaxCutoffRow(iend))
388     else
389     deallocate(gtypeMaxCutoffRow)
390     allocate(gtypeMaxCutoffRow(iend))
391     endif
392    
393    
394     #ifdef IS_MPI
395     ! We only allocate new storage if we are in MPI because Ncol /= Nrow
396 chuckv 652 if(.not.associated(groupToGtypeCol)) then
397 chuckv 651 allocate(groupToGtypeCol(jend))
398     else
399     deallocate(groupToGtypeCol)
400     allocate(groupToGtypeCol(jend))
401     end if
402    
403 tim 833 if(.not.associated(groupMaxCutoffCol)) then
404     allocate(groupMaxCutoffCol(jend))
405 chuckv 651 else
406 tim 833 deallocate(groupMaxCutoffCol)
407     allocate(groupMaxCutoffCol(jend))
408 chuckv 651 end if
409 chuckv 652 if(.not.associated(gtypeMaxCutoffCol)) then
410 chuckv 651 allocate(gtypeMaxCutoffCol(jend))
411     else
412     deallocate(gtypeMaxCutoffCol)
413     allocate(gtypeMaxCutoffCol(jend))
414     end if
415    
416     groupMaxCutoffCol = 0.0_dp
417     gtypeMaxCutoffCol = 0.0_dp
418    
419     #endif
420     groupMaxCutoffRow = 0.0_dp
421     gtypeMaxCutoffRow = 0.0_dp
422    
423    
424 gezelter 582 !! first we do a single loop over the cutoff groups to find the
425     !! largest cutoff for any atypes present in this group. We also
426     !! create gtypes at this point.
427    
428 gezelter 581 tol = 1.0d-6
429 chuckv 651 nGroupTypesRow = 0
430 tim 833 nGroupTypesCol = 0
431 gezelter 581 do i = istart, iend
432 gezelter 575 n_in_i = groupStartRow(i+1) - groupStartRow(i)
433 chuckv 651 groupMaxCutoffRow(i) = 0.0_dp
434 gezelter 581 do ia = groupStartRow(i), groupStartRow(i+1)-1
435     atom1 = groupListRow(ia)
436 gezelter 575 #ifdef IS_MPI
437 gezelter 581 me_i = atid_row(atom1)
438 gezelter 575 #else
439 gezelter 581 me_i = atid(atom1)
440     #endif
441 chuckv 651 if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
442     groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
443 gezelter 587 endif
444 gezelter 581 enddo
445 chuckv 651 if (nGroupTypesRow.eq.0) then
446     nGroupTypesRow = nGroupTypesRow + 1
447     gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
448     groupToGtypeRow(i) = nGroupTypesRow
449 gezelter 581 else
450 gezelter 587 GtypeFound = .false.
451 chuckv 651 do g = 1, nGroupTypesRow
452     if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
453     groupToGtypeRow(i) = g
454 gezelter 587 GtypeFound = .true.
455 gezelter 581 endif
456     enddo
457 gezelter 587 if (.not.GtypeFound) then
458 chuckv 651 nGroupTypesRow = nGroupTypesRow + 1
459     gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
460     groupToGtypeRow(i) = nGroupTypesRow
461 gezelter 587 endif
462 gezelter 581 endif
463 gezelter 587 enddo
464    
465 chuckv 651 #ifdef IS_MPI
466     do j = jstart, jend
467     n_in_j = groupStartCol(j+1) - groupStartCol(j)
468     groupMaxCutoffCol(j) = 0.0_dp
469     do ja = groupStartCol(j), groupStartCol(j+1)-1
470     atom1 = groupListCol(ja)
471    
472     me_j = atid_col(atom1)
473    
474     if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
475     groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
476     endif
477     enddo
478    
479     if (nGroupTypesCol.eq.0) then
480     nGroupTypesCol = nGroupTypesCol + 1
481     gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
482     groupToGtypeCol(j) = nGroupTypesCol
483     else
484     GtypeFound = .false.
485     do g = 1, nGroupTypesCol
486     if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
487     groupToGtypeCol(j) = g
488     GtypeFound = .true.
489     endif
490     enddo
491     if (.not.GtypeFound) then
492     nGroupTypesCol = nGroupTypesCol + 1
493     gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
494     groupToGtypeCol(j) = nGroupTypesCol
495     endif
496     endif
497     enddo
498    
499     #else
500     ! Set pointers to information we just found
501     nGroupTypesCol = nGroupTypesRow
502     groupToGtypeCol => groupToGtypeRow
503     gtypeMaxCutoffCol => gtypeMaxCutoffRow
504     groupMaxCutoffCol => groupMaxCutoffRow
505     #endif
506    
507 gezelter 581 !! allocate the gtypeCutoffMap here.
508 chuckv 651 allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
509 gezelter 581 !! then we do a double loop over all the group TYPES to find the cutoff
510     !! map between groups of two types
511 chuckv 651 tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
512    
513 gezelter 762 do i = 1, nGroupTypesRow
514 chuckv 651 do j = 1, nGroupTypesCol
515 gezelter 576
516 gezelter 581 select case(cutoffPolicy)
517 gezelter 582 case(TRADITIONAL_CUTOFF_POLICY)
518 chuckv 651 thisRcut = tradRcut
519 gezelter 582 case(MIX_CUTOFF_POLICY)
520 chuckv 651 thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
521 gezelter 582 case(MAX_CUTOFF_POLICY)
522 chuckv 651 thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
523 gezelter 582 case default
524     call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
525     return
526     end select
527     gtypeCutoffMap(i,j)%rcut = thisRcut
528 gezelter 762
529     if (thisRcut.gt.largestRcut) largestRcut = thisRcut
530    
531 gezelter 582 gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
532 gezelter 585
533 gezelter 762 if (.not.haveSkinThickness) then
534     skinThickness = 1.0_dp
535     endif
536    
537     gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
538    
539 chrisfen 618 ! sanity check
540    
541     if (haveDefaultCutoffs) then
542     if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
543     call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
544     endif
545     endif
546 gezelter 581 enddo
547     enddo
548 gezelter 762
549 chuckv 651 if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
550     if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
551     if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
552     #ifdef IS_MPI
553     if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
554     if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
555     #endif
556     groupMaxCutoffCol => null()
557     gtypeMaxCutoffCol => null()
558    
559 gezelter 581 haveGtypeCutoffMap = .true.
560 chrisfen 596 end subroutine createGtypeCutoffMap
561 chrisfen 578
562 gezelter 762 subroutine setCutoffs(defRcut, defRsw)
563 chrisfen 596
564 gezelter 762 real(kind=dp),intent(in) :: defRcut, defRsw
565     character(len = statusMsgSize) :: errMsg
566     integer :: localError
567    
568 chrisfen 596 defaultRcut = defRcut
569     defaultRsw = defRsw
570 gezelter 762
571     defaultDoShift = .false.
572     if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
573    
574     write(errMsg, *) &
575     'cutoffRadius and switchingRadius are set to the same', newline &
576     // tab, 'value. OOPSE will use shifted ', newline &
577     // tab, 'potentials instead of switching functions.'
578    
579     call handleInfo("setCutoffs", errMsg)
580    
581     defaultDoShift = .true.
582    
583     endif
584 chrisfen 618
585 gezelter 762 localError = 0
586     call setLJDefaultCutoff( defaultRcut, defaultDoShift )
587 gezelter 813 call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
588 gezelter 762 call setCutoffEAM( defaultRcut, localError)
589     if (localError /= 0) then
590     write(errMsg, *) 'An error has occured in setting the EAM cutoff'
591     call handleError("setCutoffs", errMsg)
592     end if
593     call set_switch(GROUP_SWITCH, defaultRsw, defaultRcut)
594 gezelter 889 call setHmatDangerousRcutValue(defaultRcut)
595 gezelter 813
596 chrisfen 618 haveDefaultCutoffs = .true.
597 gezelter 813 haveGtypeCutoffMap = .false.
598 gezelter 762 end subroutine setCutoffs
599 chrisfen 596
600 gezelter 762 subroutine cWasLame()
601    
602     VisitCutoffsAfterComputing = .true.
603     return
604    
605     end subroutine cWasLame
606    
607 chrisfen 596 subroutine setCutoffPolicy(cutPolicy)
608 gezelter 762
609 chrisfen 596 integer, intent(in) :: cutPolicy
610 gezelter 762
611 chrisfen 596 cutoffPolicy = cutPolicy
612 gezelter 762 haveCutoffPolicy = .true.
613 gezelter 813 haveGtypeCutoffMap = .false.
614 gezelter 762
615 gezelter 576 end subroutine setCutoffPolicy
616 gezelter 762
617     subroutine setElectrostaticMethod( thisESM )
618    
619     integer, intent(in) :: thisESM
620    
621     electrostaticSummationMethod = thisESM
622     haveElectrostaticSummationMethod = .true.
623 gezelter 574
624 gezelter 762 end subroutine setElectrostaticMethod
625    
626     subroutine setSkinThickness( thisSkin )
627 gezelter 574
628 gezelter 762 real(kind=dp), intent(in) :: thisSkin
629    
630     skinThickness = thisSkin
631 gezelter 813 haveSkinThickness = .true.
632     haveGtypeCutoffMap = .false.
633 gezelter 762
634     end subroutine setSkinThickness
635    
636     subroutine setSimVariables()
637     SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
638     SIM_uses_EAM = SimUsesEAM()
639     SIM_requires_postpair_calc = SimRequiresPostpairCalc()
640     SIM_requires_prepair_calc = SimRequiresPrepairCalc()
641     SIM_uses_PBC = SimUsesPBC()
642 chuckv 841 SIM_uses_SC = SimUsesSC()
643 gezelter 762
644     haveSIMvariables = .true.
645    
646     return
647     end subroutine setSimVariables
648 gezelter 117
649     subroutine doReadyCheck(error)
650     integer, intent(out) :: error
651    
652     integer :: myStatus
653    
654     error = 0
655 chrisfen 532
656 gezelter 571 if (.not. haveInteractionHash) then
657 gezelter 762 call createInteractionHash()
658 gezelter 117 endif
659    
660 gezelter 571 if (.not. haveGtypeCutoffMap) then
661 gezelter 762 call createGtypeCutoffMap()
662 gezelter 571 endif
663    
664 gezelter 762
665     if (VisitCutoffsAfterComputing) then
666     call set_switch(GROUP_SWITCH, largestRcut, largestRcut)
667 gezelter 889 call setHmatDangerousRcutValue(largestRcut)
668 gezelter 762 endif
669    
670    
671 gezelter 117 if (.not. haveSIMvariables) then
672     call setSimVariables()
673     endif
674    
675 chuckv 583 ! if (.not. haveRlist) then
676     ! write(default_error, *) 'rList has not been set in doForces!'
677     ! error = -1
678     ! return
679     ! endif
680 gezelter 117
681     if (.not. haveNeighborList) then
682     write(default_error, *) 'neighbor list has not been initialized in doForces!'
683     error = -1
684     return
685     end if
686    
687     if (.not. haveSaneForceField) then
688     write(default_error, *) 'Force Field is not sane in doForces!'
689     error = -1
690     return
691     end if
692    
693     #ifdef IS_MPI
694     if (.not. isMPISimSet()) then
695     write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
696     error = -1
697     return
698     endif
699     #endif
700     return
701     end subroutine doReadyCheck
702    
703 chrisfen 532
704 gezelter 762 subroutine init_FF(thisStat)
705 gezelter 117
706     integer, intent(out) :: thisStat
707     integer :: my_status, nMatches
708     integer, pointer :: MatchList(:) => null()
709    
710     !! assume things are copacetic, unless they aren't
711     thisStat = 0
712    
713     !! init_FF is called *after* all of the atom types have been
714     !! defined in atype_module using the new_atype subroutine.
715     !!
716     !! this will scan through the known atypes and figure out what
717     !! interactions are used by the force field.
718 chrisfen 532
719 gezelter 141 FF_uses_DirectionalAtoms = .false.
720     FF_uses_Dipoles = .false.
721     FF_uses_GayBerne = .false.
722 gezelter 117 FF_uses_EAM = .false.
723 chuckv 834 FF_uses_SC = .false.
724 chrisfen 532
725 gezelter 141 call getMatchingElementList(atypes, "is_Directional", .true., &
726     nMatches, MatchList)
727     if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
728    
729     call getMatchingElementList(atypes, "is_Dipole", .true., &
730     nMatches, MatchList)
731 gezelter 571 if (nMatches .gt. 0) FF_uses_Dipoles = .true.
732 chrisfen 523
733 gezelter 141 call getMatchingElementList(atypes, "is_GayBerne", .true., &
734     nMatches, MatchList)
735 gezelter 571 if (nMatches .gt. 0) FF_uses_GayBerne = .true.
736 chrisfen 532
737 gezelter 117 call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
738     if (nMatches .gt. 0) FF_uses_EAM = .true.
739 chrisfen 532
740 chuckv 834 call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
741     if (nMatches .gt. 0) FF_uses_SC = .true.
742 gezelter 141
743 chuckv 834
744 gezelter 117 haveSaneForceField = .true.
745 chrisfen 532
746 gezelter 117 if (FF_uses_EAM) then
747 chrisfen 532 call init_EAM_FF(my_status)
748 gezelter 117 if (my_status /= 0) then
749     write(default_error, *) "init_EAM_FF returned a bad status"
750     thisStat = -1
751     haveSaneForceField = .false.
752     return
753     end if
754     endif
755    
756     if (.not. haveNeighborList) then
757     !! Create neighbor lists
758     call expandNeighborList(nLocal, my_status)
759     if (my_Status /= 0) then
760     write(default_error,*) "SimSetup: ExpandNeighborList returned error."
761     thisStat = -1
762     return
763     endif
764     haveNeighborList = .true.
765 chrisfen 532 endif
766    
767 gezelter 117 end subroutine init_FF
768    
769 chrisfen 532
770 gezelter 117 !! Does force loop over i,j pairs. Calls do_pair to calculates forces.
771     !------------------------------------------------------------->
772 gezelter 246 subroutine do_force_loop(q, q_group, A, eFrame, f, t, tau, pot, &
773 gezelter 117 do_pot_c, do_stress_c, error)
774     !! Position array provided by C, dimensioned by getNlocal
775     real ( kind = dp ), dimension(3, nLocal) :: q
776     !! molecular center-of-mass position array
777     real ( kind = dp ), dimension(3, nGroups) :: q_group
778     !! Rotation Matrix for each long range particle in simulation.
779     real( kind = dp), dimension(9, nLocal) :: A
780     !! Unit vectors for dipoles (lab frame)
781 gezelter 246 real( kind = dp ), dimension(9,nLocal) :: eFrame
782 gezelter 117 !! Force array provided by C, dimensioned by getNlocal
783     real ( kind = dp ), dimension(3,nLocal) :: f
784     !! Torsion array provided by C, dimensioned by getNlocal
785     real( kind = dp ), dimension(3,nLocal) :: t
786    
787     !! Stress Tensor
788     real( kind = dp), dimension(9) :: tau
789 gezelter 662 real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
790 gezelter 117 logical ( kind = 2) :: do_pot_c, do_stress_c
791     logical :: do_pot
792     logical :: do_stress
793     logical :: in_switching_region
794     #ifdef IS_MPI
795 gezelter 662 real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
796 gezelter 117 integer :: nAtomsInRow
797     integer :: nAtomsInCol
798     integer :: nprocs
799     integer :: nGroupsInRow
800     integer :: nGroupsInCol
801     #endif
802     integer :: natoms
803     logical :: update_nlist
804     integer :: i, j, jstart, jend, jnab
805     integer :: istart, iend
806     integer :: ia, jb, atom1, atom2
807     integer :: nlist
808     real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
809     real( kind = DP ) :: sw, dswdr, swderiv, mf
810 chrisfen 699 real( kind = DP ) :: rVal
811 gezelter 117 real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
812     real(kind=dp) :: rfpot, mu_i, virial
813 gezelter 762 real(kind=dp):: rCut
814 gezelter 117 integer :: me_i, me_j, n_in_i, n_in_j
815     logical :: is_dp_i
816     integer :: neighborListSize
817     integer :: listerror, error
818     integer :: localError
819     integer :: propPack_i, propPack_j
820     integer :: loopStart, loopEnd, loop
821 gezelter 571 integer :: iHash
822 chrisfen 699 integer :: i1
823 chuckv 624
824 chrisfen 532
825 gezelter 117 !! initialize local variables
826 chrisfen 532
827 gezelter 117 #ifdef IS_MPI
828     pot_local = 0.0_dp
829     nAtomsInRow = getNatomsInRow(plan_atom_row)
830     nAtomsInCol = getNatomsInCol(plan_atom_col)
831     nGroupsInRow = getNgroupsInRow(plan_group_row)
832     nGroupsInCol = getNgroupsInCol(plan_group_col)
833     #else
834     natoms = nlocal
835     #endif
836 chrisfen 532
837 gezelter 117 call doReadyCheck(localError)
838     if ( localError .ne. 0 ) then
839     call handleError("do_force_loop", "Not Initialized")
840     error = -1
841     return
842     end if
843     call zero_work_arrays()
844 chrisfen 532
845 gezelter 117 do_pot = do_pot_c
846     do_stress = do_stress_c
847 chrisfen 532
848 gezelter 117 ! Gather all information needed by all force loops:
849 chrisfen 532
850 gezelter 117 #ifdef IS_MPI
851 chrisfen 532
852 gezelter 117 call gather(q, q_Row, plan_atom_row_3d)
853     call gather(q, q_Col, plan_atom_col_3d)
854    
855     call gather(q_group, q_group_Row, plan_group_row_3d)
856     call gather(q_group, q_group_Col, plan_group_col_3d)
857 chrisfen 532
858 gezelter 141 if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
859 gezelter 246 call gather(eFrame, eFrame_Row, plan_atom_row_rotation)
860     call gather(eFrame, eFrame_Col, plan_atom_col_rotation)
861 chrisfen 532
862 gezelter 117 call gather(A, A_Row, plan_atom_row_rotation)
863     call gather(A, A_Col, plan_atom_col_rotation)
864     endif
865 chrisfen 532
866 gezelter 117 #endif
867 chrisfen 532
868 gezelter 117 !! Begin force loop timing:
869     #ifdef PROFILE
870     call cpu_time(forceTimeInitial)
871     nloops = nloops + 1
872     #endif
873 chrisfen 532
874 gezelter 117 loopEnd = PAIR_LOOP
875     if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
876     loopStart = PREPAIR_LOOP
877     else
878     loopStart = PAIR_LOOP
879     endif
880    
881     do loop = loopStart, loopEnd
882    
883     ! See if we need to update neighbor lists
884     ! (but only on the first time through):
885     if (loop .eq. loopStart) then
886     #ifdef IS_MPI
887 gezelter 762 call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
888 chrisfen 532 update_nlist)
889 gezelter 117 #else
890 gezelter 762 call checkNeighborList(nGroups, q_group, skinThickness, &
891 chrisfen 532 update_nlist)
892 gezelter 117 #endif
893     endif
894 chrisfen 532
895 gezelter 117 if (update_nlist) then
896     !! save current configuration and construct neighbor list
897     #ifdef IS_MPI
898     call saveNeighborList(nGroupsInRow, q_group_row)
899     #else
900     call saveNeighborList(nGroups, q_group)
901     #endif
902     neighborListSize = size(list)
903     nlist = 0
904     endif
905 chrisfen 532
906 gezelter 117 istart = 1
907     #ifdef IS_MPI
908     iend = nGroupsInRow
909     #else
910     iend = nGroups - 1
911     #endif
912     outer: do i = istart, iend
913    
914     if (update_nlist) point(i) = nlist + 1
915 chrisfen 532
916 gezelter 117 n_in_i = groupStartRow(i+1) - groupStartRow(i)
917 chrisfen 532
918 gezelter 117 if (update_nlist) then
919     #ifdef IS_MPI
920     jstart = 1
921     jend = nGroupsInCol
922     #else
923     jstart = i+1
924     jend = nGroups
925     #endif
926     else
927     jstart = point(i)
928     jend = point(i+1) - 1
929     ! make sure group i has neighbors
930     if (jstart .gt. jend) cycle outer
931     endif
932 chrisfen 532
933 gezelter 117 do jnab = jstart, jend
934     if (update_nlist) then
935     j = jnab
936     else
937     j = list(jnab)
938     endif
939    
940     #ifdef IS_MPI
941 chuckv 567 me_j = atid_col(j)
942 gezelter 117 call get_interatomic_vector(q_group_Row(:,i), &
943     q_group_Col(:,j), d_grp, rgrpsq)
944     #else
945 chuckv 567 me_j = atid(j)
946 gezelter 117 call get_interatomic_vector(q_group(:,i), &
947     q_group(:,j), d_grp, rgrpsq)
948 chrisfen 618 #endif
949 gezelter 117
950 chuckv 651 if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
951 gezelter 117 if (update_nlist) then
952     nlist = nlist + 1
953 chrisfen 532
954 gezelter 117 if (nlist > neighborListSize) then
955     #ifdef IS_MPI
956     call expandNeighborList(nGroupsInRow, listerror)
957     #else
958     call expandNeighborList(nGroups, listerror)
959     #endif
960     if (listerror /= 0) then
961     error = -1
962     write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
963     return
964     end if
965     neighborListSize = size(list)
966     endif
967 chrisfen 532
968 gezelter 117 list(nlist) = j
969     endif
970 gezelter 762
971    
972 chrisfen 708
973     if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
974 chrisfen 532
975 gezelter 762 rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
976 chrisfen 708 if (loop .eq. PAIR_LOOP) then
977     vij = 0.0d0
978     fij(1:3) = 0.0d0
979     endif
980    
981 gezelter 762 call get_switch(rgrpsq, sw, dswdr, rgrp, &
982     group_switch, in_switching_region)
983 chrisfen 708
984     n_in_j = groupStartCol(j+1) - groupStartCol(j)
985    
986     do ia = groupStartRow(i), groupStartRow(i+1)-1
987 chrisfen 703
988 chrisfen 708 atom1 = groupListRow(ia)
989    
990     inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
991    
992     atom2 = groupListCol(jb)
993    
994     if (skipThisPair(atom1, atom2)) cycle inner
995    
996     if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
997     d_atm(1:3) = d_grp(1:3)
998     ratmsq = rgrpsq
999     else
1000 gezelter 117 #ifdef IS_MPI
1001 chrisfen 708 call get_interatomic_vector(q_Row(:,atom1), &
1002     q_Col(:,atom2), d_atm, ratmsq)
1003 gezelter 117 #else
1004 chrisfen 708 call get_interatomic_vector(q(:,atom1), &
1005     q(:,atom2), d_atm, ratmsq)
1006 gezelter 117 #endif
1007 chrisfen 708 endif
1008    
1009     if (loop .eq. PREPAIR_LOOP) then
1010 gezelter 117 #ifdef IS_MPI
1011 chrisfen 708 call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1012 gezelter 762 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1013 chrisfen 708 eFrame, A, f, t, pot_local)
1014 gezelter 117 #else
1015 chrisfen 708 call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1016 gezelter 762 rgrpsq, d_grp, rCut, do_pot, do_stress, &
1017 chrisfen 708 eFrame, A, f, t, pot)
1018 gezelter 117 #endif
1019 chrisfen 708 else
1020 gezelter 117 #ifdef IS_MPI
1021 chrisfen 708 call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1022     do_pot, eFrame, A, f, t, pot_local, vpair, &
1023 gezelter 762 fpair, d_grp, rgrp, rCut)
1024 gezelter 117 #else
1025 chrisfen 708 call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1026     do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1027 gezelter 762 d_grp, rgrp, rCut)
1028 gezelter 117 #endif
1029 chrisfen 708 vij = vij + vpair
1030     fij(1:3) = fij(1:3) + fpair(1:3)
1031     endif
1032     enddo inner
1033     enddo
1034 gezelter 117
1035 chrisfen 708 if (loop .eq. PAIR_LOOP) then
1036     if (in_switching_region) then
1037     swderiv = vij*dswdr/rgrp
1038     fij(1) = fij(1) + swderiv*d_grp(1)
1039     fij(2) = fij(2) + swderiv*d_grp(2)
1040     fij(3) = fij(3) + swderiv*d_grp(3)
1041    
1042     do ia=groupStartRow(i), groupStartRow(i+1)-1
1043     atom1=groupListRow(ia)
1044     mf = mfactRow(atom1)
1045 gezelter 117 #ifdef IS_MPI
1046 chrisfen 708 f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1047     f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1048     f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1049 gezelter 117 #else
1050 chrisfen 708 f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1051     f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1052     f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1053 gezelter 117 #endif
1054 chrisfen 708 enddo
1055    
1056     do jb=groupStartCol(j), groupStartCol(j+1)-1
1057     atom2=groupListCol(jb)
1058     mf = mfactCol(atom2)
1059 gezelter 117 #ifdef IS_MPI
1060 chrisfen 708 f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1061     f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1062     f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1063 gezelter 117 #else
1064 chrisfen 708 f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1065     f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1066     f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1067 gezelter 117 #endif
1068 chrisfen 708 enddo
1069     endif
1070    
1071     if (do_stress) call add_stress_tensor(d_grp, fij)
1072 gezelter 117 endif
1073     endif
1074 chrisfen 708 endif
1075 gezelter 117 enddo
1076 chrisfen 708
1077 gezelter 117 enddo outer
1078 chrisfen 532
1079 gezelter 117 if (update_nlist) then
1080     #ifdef IS_MPI
1081     point(nGroupsInRow + 1) = nlist + 1
1082     #else
1083     point(nGroups) = nlist + 1
1084     #endif
1085     if (loop .eq. PREPAIR_LOOP) then
1086     ! we just did the neighbor list update on the first
1087     ! pass, so we don't need to do it
1088     ! again on the second pass
1089     update_nlist = .false.
1090     endif
1091     endif
1092 chrisfen 532
1093 gezelter 117 if (loop .eq. PREPAIR_LOOP) then
1094     call do_preforce(nlocal, pot)
1095     endif
1096 chrisfen 532
1097 gezelter 117 enddo
1098 chrisfen 532
1099 gezelter 117 !! Do timing
1100     #ifdef PROFILE
1101     call cpu_time(forceTimeFinal)
1102     forceTime = forceTime + forceTimeFinal - forceTimeInitial
1103     #endif
1104 chrisfen 532
1105 gezelter 117 #ifdef IS_MPI
1106     !!distribute forces
1107 chrisfen 532
1108 gezelter 117 f_temp = 0.0_dp
1109     call scatter(f_Row,f_temp,plan_atom_row_3d)
1110     do i = 1,nlocal
1111     f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
1112     end do
1113 chrisfen 532
1114 gezelter 117 f_temp = 0.0_dp
1115     call scatter(f_Col,f_temp,plan_atom_col_3d)
1116     do i = 1,nlocal
1117     f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
1118     end do
1119 chrisfen 532
1120 gezelter 141 if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
1121 gezelter 117 t_temp = 0.0_dp
1122     call scatter(t_Row,t_temp,plan_atom_row_3d)
1123     do i = 1,nlocal
1124     t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
1125     end do
1126     t_temp = 0.0_dp
1127     call scatter(t_Col,t_temp,plan_atom_col_3d)
1128 chrisfen 532
1129 gezelter 117 do i = 1,nlocal
1130     t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
1131     end do
1132     endif
1133 chrisfen 532
1134 gezelter 117 if (do_pot) then
1135     ! scatter/gather pot_row into the members of my column
1136 gezelter 662 do i = 1,LR_POT_TYPES
1137 chuckv 657 call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1138     end do
1139 gezelter 117 ! scatter/gather pot_local into all other procs
1140     ! add resultant to get total pot
1141     do i = 1, nlocal
1142 gezelter 662 pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1143     + pot_Temp(1:LR_POT_TYPES,i)
1144 gezelter 117 enddo
1145 chrisfen 532
1146 gezelter 117 pot_Temp = 0.0_DP
1147 gezelter 662 do i = 1,LR_POT_TYPES
1148 chuckv 657 call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1149     end do
1150 gezelter 117 do i = 1, nlocal
1151 gezelter 662 pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1152     + pot_Temp(1:LR_POT_TYPES,i)
1153 gezelter 117 enddo
1154 chrisfen 532
1155 gezelter 117 endif
1156     #endif
1157 chrisfen 532
1158 chrisfen 691 if (SIM_requires_postpair_calc) then
1159 chrisfen 695 do i = 1, nlocal
1160    
1161     ! we loop only over the local atoms, so we don't need row and column
1162     ! lookups for the types
1163 chrisfen 699
1164 chrisfen 691 me_i = atid(i)
1165    
1166 chrisfen 695 ! is the atom electrostatic? See if it would have an
1167     ! electrostatic interaction with itself
1168     iHash = InteractionHash(me_i,me_i)
1169 chrisfen 699
1170 chrisfen 691 if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1171 gezelter 117 #ifdef IS_MPI
1172 chrisfen 703 call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1173 chrisfen 695 t, do_pot)
1174 gezelter 117 #else
1175 chrisfen 703 call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1176 chrisfen 695 t, do_pot)
1177 gezelter 117 #endif
1178 chrisfen 691 endif
1179 chrisfen 699
1180 chrisfen 703
1181 chrisfen 708 if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1182 chrisfen 699
1183 chrisfen 703 ! loop over the excludes to accumulate RF stuff we've
1184     ! left out of the normal pair loop
1185    
1186     do i1 = 1, nSkipsForAtom(i)
1187     j = skipsForAtom(i, i1)
1188    
1189     ! prevent overcounting of the skips
1190     if (i.lt.j) then
1191     call get_interatomic_vector(q(:,i), &
1192     q(:,j), d_atm, ratmsq)
1193     rVal = dsqrt(ratmsq)
1194     call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1195     in_switching_region)
1196 chrisfen 699 #ifdef IS_MPI
1197 chrisfen 703 call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1198     vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1199 chrisfen 699 #else
1200 chrisfen 703 call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1201     vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1202 chrisfen 699 #endif
1203 chrisfen 703 endif
1204     enddo
1205 chrisfen 708 endif
1206 chrisfen 703 enddo
1207 gezelter 117 endif
1208 chrisfen 691
1209 gezelter 117 #ifdef IS_MPI
1210 chrisfen 695
1211 gezelter 117 if (do_pot) then
1212 chrisfen 695 call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1213     mpi_comm_world,mpi_err)
1214 gezelter 117 endif
1215 chrisfen 695
1216 gezelter 117 if (do_stress) then
1217     call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1218     mpi_comm_world,mpi_err)
1219     call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1220     mpi_comm_world,mpi_err)
1221     endif
1222 chrisfen 695
1223 gezelter 117 #else
1224 chrisfen 695
1225 gezelter 117 if (do_stress) then
1226     tau = tau_Temp
1227     virial = virial_Temp
1228     endif
1229 chrisfen 695
1230 gezelter 117 #endif
1231 chrisfen 695
1232 gezelter 117 end subroutine do_force_loop
1233 chrisfen 532
1234 gezelter 117 subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
1235 gezelter 762 eFrame, A, f, t, pot, vpair, fpair, d_grp, r_grp, rCut)
1236 gezelter 117
1237 chuckv 656 real( kind = dp ) :: vpair, sw
1238 gezelter 662 real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1239 gezelter 117 real( kind = dp ), dimension(3) :: fpair
1240     real( kind = dp ), dimension(nLocal) :: mfact
1241 gezelter 246 real( kind = dp ), dimension(9,nLocal) :: eFrame
1242 gezelter 117 real( kind = dp ), dimension(9,nLocal) :: A
1243     real( kind = dp ), dimension(3,nLocal) :: f
1244     real( kind = dp ), dimension(3,nLocal) :: t
1245    
1246     logical, intent(inout) :: do_pot
1247     integer, intent(in) :: i, j
1248     real ( kind = dp ), intent(inout) :: rijsq
1249 chrisfen 695 real ( kind = dp ), intent(inout) :: r_grp
1250 gezelter 117 real ( kind = dp ), intent(inout) :: d(3)
1251 chrisfen 695 real ( kind = dp ), intent(inout) :: d_grp(3)
1252 gezelter 762 real ( kind = dp ), intent(inout) :: rCut
1253 chrisfen 695 real ( kind = dp ) :: r
1254 gezelter 117 integer :: me_i, me_j
1255    
1256 gezelter 571 integer :: iHash
1257 gezelter 560
1258 gezelter 117 r = sqrt(rijsq)
1259     vpair = 0.0d0
1260     fpair(1:3) = 0.0d0
1261    
1262     #ifdef IS_MPI
1263     me_i = atid_row(i)
1264     me_j = atid_col(j)
1265     #else
1266     me_i = atid(i)
1267     me_j = atid(j)
1268     #endif
1269 gezelter 202
1270 gezelter 571 iHash = InteractionHash(me_i, me_j)
1271 chrisfen 703
1272     if ( iand(iHash, LJ_PAIR).ne.0 ) then
1273 gezelter 762 call do_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1274 chrisfen 703 pot(VDW_POT), f, do_pot)
1275 gezelter 117 endif
1276 chrisfen 532
1277 chrisfen 703 if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1278 gezelter 762 call doElectrostaticPair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1279 chrisfen 712 pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1280 chrisfen 703 endif
1281    
1282     if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1283     call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1284     pot(HB_POT), A, f, t, do_pot)
1285     endif
1286    
1287     if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1288     call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1289     pot(HB_POT), A, f, t, do_pot)
1290     endif
1291    
1292     if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1293     call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1294     pot(VDW_POT), A, f, t, do_pot)
1295     endif
1296    
1297     if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1298 gezelter 762 call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1299 chrisfen 703 pot(VDW_POT), A, f, t, do_pot)
1300     endif
1301    
1302     if ( iand(iHash, EAM_PAIR).ne.0 ) then
1303     call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1304     pot(METALLIC_POT), f, do_pot)
1305     endif
1306    
1307     if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1308     call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1309     pot(VDW_POT), A, f, t, do_pot)
1310     endif
1311    
1312     if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1313     call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1314     pot(VDW_POT), A, f, t, do_pot)
1315     endif
1316 chuckv 733
1317     if ( iand(iHash, SC_PAIR).ne.0 ) then
1318 gezelter 762 call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1319 chuckv 733 pot(METALLIC_POT), f, do_pot)
1320     endif
1321    
1322    
1323 chrisfen 703
1324 gezelter 117 end subroutine do_pair
1325    
1326 gezelter 762 subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1327 gezelter 246 do_pot, do_stress, eFrame, A, f, t, pot)
1328 gezelter 117
1329 chuckv 656 real( kind = dp ) :: sw
1330 gezelter 662 real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1331 chrisfen 532 real( kind = dp ), dimension(9,nLocal) :: eFrame
1332     real (kind=dp), dimension(9,nLocal) :: A
1333     real (kind=dp), dimension(3,nLocal) :: f
1334     real (kind=dp), dimension(3,nLocal) :: t
1335 gezelter 117
1336 chrisfen 532 logical, intent(inout) :: do_pot, do_stress
1337     integer, intent(in) :: i, j
1338 gezelter 762 real ( kind = dp ), intent(inout) :: rijsq, rcijsq, rCut
1339 chrisfen 532 real ( kind = dp ) :: r, rc
1340     real ( kind = dp ), intent(inout) :: d(3), dc(3)
1341    
1342 gezelter 571 integer :: me_i, me_j, iHash
1343 chrisfen 532
1344 gezelter 591 r = sqrt(rijsq)
1345    
1346 gezelter 117 #ifdef IS_MPI
1347 chrisfen 532 me_i = atid_row(i)
1348     me_j = atid_col(j)
1349 gezelter 117 #else
1350 chrisfen 532 me_i = atid(i)
1351     me_j = atid(j)
1352 gezelter 117 #endif
1353 chrisfen 532
1354 gezelter 571 iHash = InteractionHash(me_i, me_j)
1355 chrisfen 532
1356 gezelter 571 if ( iand(iHash, EAM_PAIR).ne.0 ) then
1357 gezelter 762 call calc_EAM_prepair_rho(i, j, d, r, rijsq)
1358 chrisfen 532 endif
1359 chuckv 733
1360     if ( iand(iHash, SC_PAIR).ne.0 ) then
1361 gezelter 762 call calc_SC_prepair_rho(i, j, d, r, rijsq, rcut )
1362 chuckv 733 endif
1363 gezelter 560
1364 chrisfen 532 end subroutine do_prepair
1365    
1366    
1367     subroutine do_preforce(nlocal,pot)
1368     integer :: nlocal
1369 gezelter 662 real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1370 chrisfen 532
1371     if (FF_uses_EAM .and. SIM_uses_EAM) then
1372 gezelter 662 call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1373 chrisfen 532 endif
1374 chuckv 733 if (FF_uses_SC .and. SIM_uses_SC) then
1375     call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1376     endif
1377 chrisfen 532
1378    
1379     end subroutine do_preforce
1380    
1381    
1382     subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1383    
1384     real (kind = dp), dimension(3) :: q_i
1385     real (kind = dp), dimension(3) :: q_j
1386     real ( kind = dp ), intent(out) :: r_sq
1387     real( kind = dp ) :: d(3), scaled(3)
1388     integer i
1389    
1390     d(1:3) = q_j(1:3) - q_i(1:3)
1391    
1392     ! Wrap back into periodic box if necessary
1393     if ( SIM_uses_PBC ) then
1394    
1395     if( .not.boxIsOrthorhombic ) then
1396     ! calc the scaled coordinates.
1397    
1398     scaled = matmul(HmatInv, d)
1399    
1400     ! wrap the scaled coordinates
1401    
1402     scaled = scaled - anint(scaled)
1403    
1404    
1405     ! calc the wrapped real coordinates from the wrapped scaled
1406     ! coordinates
1407    
1408     d = matmul(Hmat,scaled)
1409    
1410     else
1411     ! calc the scaled coordinates.
1412    
1413     do i = 1, 3
1414     scaled(i) = d(i) * HmatInv(i,i)
1415    
1416     ! wrap the scaled coordinates
1417    
1418     scaled(i) = scaled(i) - anint(scaled(i))
1419    
1420     ! calc the wrapped real coordinates from the wrapped scaled
1421     ! coordinates
1422    
1423     d(i) = scaled(i)*Hmat(i,i)
1424     enddo
1425     endif
1426    
1427     endif
1428    
1429     r_sq = dot_product(d,d)
1430    
1431     end subroutine get_interatomic_vector
1432    
1433     subroutine zero_work_arrays()
1434    
1435 gezelter 117 #ifdef IS_MPI
1436    
1437 chrisfen 532 q_Row = 0.0_dp
1438     q_Col = 0.0_dp
1439    
1440     q_group_Row = 0.0_dp
1441     q_group_Col = 0.0_dp
1442    
1443     eFrame_Row = 0.0_dp
1444     eFrame_Col = 0.0_dp
1445    
1446     A_Row = 0.0_dp
1447     A_Col = 0.0_dp
1448    
1449     f_Row = 0.0_dp
1450     f_Col = 0.0_dp
1451     f_Temp = 0.0_dp
1452    
1453     t_Row = 0.0_dp
1454     t_Col = 0.0_dp
1455     t_Temp = 0.0_dp
1456    
1457     pot_Row = 0.0_dp
1458     pot_Col = 0.0_dp
1459     pot_Temp = 0.0_dp
1460    
1461 gezelter 117 #endif
1462 chrisfen 532
1463     if (FF_uses_EAM .and. SIM_uses_EAM) then
1464     call clean_EAM()
1465     endif
1466    
1467     tau_Temp = 0.0_dp
1468     virial_Temp = 0.0_dp
1469     end subroutine zero_work_arrays
1470    
1471     function skipThisPair(atom1, atom2) result(skip_it)
1472     integer, intent(in) :: atom1
1473     integer, intent(in), optional :: atom2
1474     logical :: skip_it
1475     integer :: unique_id_1, unique_id_2
1476     integer :: me_i,me_j
1477     integer :: i
1478    
1479     skip_it = .false.
1480    
1481     !! there are a number of reasons to skip a pair or a particle
1482     !! mostly we do this to exclude atoms who are involved in short
1483     !! range interactions (bonds, bends, torsions), but we also need
1484     !! to exclude some overcounted interactions that result from
1485     !! the parallel decomposition
1486    
1487 gezelter 117 #ifdef IS_MPI
1488 chrisfen 532 !! in MPI, we have to look up the unique IDs for each atom
1489     unique_id_1 = AtomRowToGlobal(atom1)
1490 gezelter 117 #else
1491 chrisfen 532 !! in the normal loop, the atom numbers are unique
1492     unique_id_1 = atom1
1493 gezelter 117 #endif
1494 chrisfen 532
1495     !! We were called with only one atom, so just check the global exclude
1496     !! list for this atom
1497     if (.not. present(atom2)) then
1498     do i = 1, nExcludes_global
1499     if (excludesGlobal(i) == unique_id_1) then
1500     skip_it = .true.
1501     return
1502     end if
1503     end do
1504     return
1505     end if
1506    
1507 gezelter 117 #ifdef IS_MPI
1508 chrisfen 532 unique_id_2 = AtomColToGlobal(atom2)
1509 gezelter 117 #else
1510 chrisfen 532 unique_id_2 = atom2
1511 gezelter 117 #endif
1512 chrisfen 532
1513 gezelter 117 #ifdef IS_MPI
1514 chrisfen 532 !! this situation should only arise in MPI simulations
1515     if (unique_id_1 == unique_id_2) then
1516     skip_it = .true.
1517     return
1518     end if
1519    
1520     !! this prevents us from doing the pair on multiple processors
1521     if (unique_id_1 < unique_id_2) then
1522     if (mod(unique_id_1 + unique_id_2,2) == 0) then
1523     skip_it = .true.
1524     return
1525     endif
1526     else
1527     if (mod(unique_id_1 + unique_id_2,2) == 1) then
1528     skip_it = .true.
1529     return
1530     endif
1531     endif
1532 gezelter 117 #endif
1533 chrisfen 532
1534     !! the rest of these situations can happen in all simulations:
1535     do i = 1, nExcludes_global
1536     if ((excludesGlobal(i) == unique_id_1) .or. &
1537     (excludesGlobal(i) == unique_id_2)) then
1538     skip_it = .true.
1539     return
1540     endif
1541     enddo
1542    
1543     do i = 1, nSkipsForAtom(atom1)
1544     if (skipsForAtom(atom1, i) .eq. unique_id_2) then
1545     skip_it = .true.
1546     return
1547     endif
1548     end do
1549    
1550     return
1551     end function skipThisPair
1552    
1553     function FF_UsesDirectionalAtoms() result(doesit)
1554     logical :: doesit
1555 gezelter 571 doesit = FF_uses_DirectionalAtoms
1556 chrisfen 532 end function FF_UsesDirectionalAtoms
1557    
1558     function FF_RequiresPrepairCalc() result(doesit)
1559     logical :: doesit
1560 chuckv 733 doesit = FF_uses_EAM .or. FF_uses_SC &
1561     .or. FF_uses_MEAM
1562 chrisfen 532 end function FF_RequiresPrepairCalc
1563    
1564 gezelter 117 #ifdef PROFILE
1565 chrisfen 532 function getforcetime() result(totalforcetime)
1566     real(kind=dp) :: totalforcetime
1567     totalforcetime = forcetime
1568     end function getforcetime
1569 gezelter 117 #endif
1570    
1571 chrisfen 532 !! This cleans componets of force arrays belonging only to fortran
1572    
1573     subroutine add_stress_tensor(dpair, fpair)
1574    
1575     real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1576    
1577     ! because the d vector is the rj - ri vector, and
1578     ! because fx, fy, fz are the force on atom i, we need a
1579     ! negative sign here:
1580    
1581     tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1582     tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1583     tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1584     tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1585     tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1586     tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1587     tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1588     tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1589     tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1590    
1591     virial_Temp = virial_Temp + &
1592     (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1593    
1594     end subroutine add_stress_tensor
1595    
1596 gezelter 117 end module doForces