ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/UseTheForce/doForces.F90
Revision: 1596
Committed: Mon Jul 25 17:30:53 2011 UTC (14 years ago) by gezelter
File size: 68259 byte(s)
Log Message:
Updated the BlockSnapshotManager to use a specified memory footprint
in constructor and not to rely on physmem and residentMem to figure
out free memory. DynamicProps is the only program that uses the
BlockSnapshotManager, so substantial changes were needed there as
well.


File Contents

# User Rev Content
1 gezelter 246 !!
2 chuckv 1388 !! Copyright (c) 2005, 2009 The University of Notre Dame. All Rights Reserved.
3 gezelter 246 !!
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 gezelter 1390 !! 1. Redistributions of source code must retain the above copyright
10 gezelter 246 !! notice, this list of conditions and the following disclaimer.
11     !!
12 gezelter 1390 !! 2. Redistributions in binary form must reproduce the above copyright
13 gezelter 246 !! notice, this list of conditions and the following disclaimer in the
14     !! documentation and/or other materials provided with the
15     !! distribution.
16     !!
17     !! This software is provided "AS IS," without a warranty of any
18     !! kind. All express or implied conditions, representations and
19     !! warranties, including any implied warranty of merchantability,
20     !! fitness for a particular purpose or non-infringement, are hereby
21     !! excluded. The University of Notre Dame and its licensors shall not
22     !! be liable for any damages suffered by licensee as a result of
23     !! using, modifying or distributing the software or its
24     !! derivatives. In no event will the University of Notre Dame or its
25     !! licensors be liable for any lost revenue, profit or data, or for
26     !! direct, indirect, special, consequential, incidental or punitive
27     !! damages, however caused and regardless of the theory of liability,
28     !! arising out of the use of or inability to use software, even if the
29     !! University of Notre Dame has been advised of the possibility of
30     !! such damages.
31     !!
32 gezelter 1390 !! SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33     !! research, please cite the appropriate papers when you publish your
34     !! work. Good starting points are:
35     !!
36     !! [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37     !! [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38     !! [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39     !! [4] Vardeman & Gezelter, in progress (2009).
40     !!
41 gezelter 246
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 gezelter 1442 !! @version $Id$, $Date$, $Name: not supported by cvs2svn $, $Revision$
49 gezelter 117
50 gezelter 246
51 gezelter 117 module doForces
52     use force_globals
53 gezelter 1286 use fForceOptions
54 gezelter 117 use simulation
55     use definitions
56     use atype_module
57     use switcheroo
58     use neighborLists
59     use lj
60 gezelter 246 use sticky
61 gezelter 401 use electrostatic_module
62 gezelter 676 use gayberne
63 chrisfen 143 use shapes
64 gezelter 117 use vector_class
65     use eam
66 chuckv 1162 use MetalNonMetal
67 chuckv 733 use suttonchen
68 gezelter 117 use status
69     #ifdef IS_MPI
70     use mpiSimulation
71     #endif
72    
73     implicit none
74     PRIVATE
75    
76     #define __FORTRAN90
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     INTEGER, PARAMETER:: PREPAIR_LOOP = 1
82     INTEGER, PARAMETER:: PAIR_LOOP = 2
83    
84     logical, save :: haveNeighborList = .false.
85     logical, save :: haveSIMvariables = .false.
86     logical, save :: haveSaneForceField = .false.
87 gezelter 571 logical, save :: haveInteractionHash = .false.
88     logical, save :: haveGtypeCutoffMap = .false.
89 chrisfen 618 logical, save :: haveDefaultCutoffs = .false.
90 gezelter 762 logical, save :: haveSkinThickness = .false.
91     logical, save :: haveElectrostaticSummationMethod = .false.
92     logical, save :: haveCutoffPolicy = .false.
93     logical, save :: VisitCutoffsAfterComputing = .false.
94 chrisfen 998 logical, save :: do_box_dipole = .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 chuckv 1162 logical, save :: FF_uses_MNM
102 chuckv 733
103 gezelter 141
104     logical, save :: SIM_uses_DirectionalAtoms
105     logical, save :: SIM_uses_EAM
106 chuckv 733 logical, save :: SIM_uses_SC
107 chuckv 1162 logical, save :: SIM_uses_MNM
108 gezelter 117 logical, save :: SIM_requires_postpair_calc
109     logical, save :: SIM_requires_prepair_calc
110     logical, save :: SIM_uses_PBC
111 gezelter 1126 logical, save :: SIM_uses_AtomicVirial
112 gezelter 117
113 chrisfen 607 integer, save :: electrostaticSummationMethod
114 gezelter 762 integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
115 chrisfen 580
116 gezelter 762 real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut
117     real(kind=dp), save :: skinThickness
118 chrisfen 1129 logical, save :: defaultDoShiftPot
119     logical, save :: defaultDoShiftFrc
120 gezelter 762
121 gezelter 117 public :: init_FF
122 gezelter 762 public :: setCutoffs
123     public :: cWasLame
124     public :: setElectrostaticMethod
125 chrisfen 998 public :: setBoxDipole
126     public :: getBoxDipole
127 gezelter 762 public :: setCutoffPolicy
128     public :: setSkinThickness
129 gezelter 117 public :: do_force_loop
130    
131     #ifdef PROFILE
132     public :: getforcetime
133     real, save :: forceTime = 0
134     real :: forceTimeInitial, forceTimeFinal
135     integer :: nLoops
136     #endif
137 chuckv 561
138 gezelter 571 !! Variables for cutoff mapping and interaction mapping
139     ! Bit hash to determine pair-pair interactions.
140     integer, dimension(:,:), allocatable :: InteractionHash
141     real(kind=dp), dimension(:), allocatable :: atypeMaxCutoff
142 chuckv 651 real(kind=dp), dimension(:), allocatable, target :: groupMaxCutoffRow
143     real(kind=dp), dimension(:), pointer :: groupMaxCutoffCol
144    
145     integer, dimension(:), allocatable, target :: groupToGtypeRow
146     integer, dimension(:), pointer :: groupToGtypeCol => null()
147    
148     real(kind=dp), dimension(:), allocatable,target :: gtypeMaxCutoffRow
149     real(kind=dp), dimension(:), pointer :: gtypeMaxCutoffCol
150 gezelter 571 type ::gtypeCutoffs
151     real(kind=dp) :: rcut
152     real(kind=dp) :: rcutsq
153     real(kind=dp) :: rlistsq
154     end type gtypeCutoffs
155     type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
156 gezelter 574
157 chrisfen 998 real(kind=dp), dimension(3) :: boxDipole
158 gezelter 939
159 gezelter 117 contains
160    
161 gezelter 762 subroutine createInteractionHash()
162 chuckv 561 integer :: nAtypes
163     integer :: i
164     integer :: j
165 gezelter 571 integer :: iHash
166 tim 568 !! Test Types
167 chuckv 561 logical :: i_is_LJ
168     logical :: i_is_Elect
169     logical :: i_is_Sticky
170     logical :: i_is_StickyP
171     logical :: i_is_GB
172     logical :: i_is_EAM
173     logical :: i_is_Shape
174 chuckv 733 logical :: i_is_SC
175 chuckv 561 logical :: j_is_LJ
176     logical :: j_is_Elect
177     logical :: j_is_Sticky
178     logical :: j_is_StickyP
179     logical :: j_is_GB
180     logical :: j_is_EAM
181     logical :: j_is_Shape
182 chuckv 733 logical :: j_is_SC
183 gezelter 576 real(kind=dp) :: myRcut
184    
185 chuckv 561 if (.not. associated(atypes)) then
186 gezelter 762 call handleError("doForces", "atypes was not present before call of createInteractionHash!")
187 chuckv 561 return
188     endif
189    
190     nAtypes = getSize(atypes)
191    
192     if (nAtypes == 0) then
193 gezelter 762 call handleError("doForces", "nAtypes was zero during call of createInteractionHash!")
194 chuckv 561 return
195     end if
196 chrisfen 532
197 chuckv 570 if (.not. allocated(InteractionHash)) then
198     allocate(InteractionHash(nAtypes,nAtypes))
199 chuckv 655 else
200     deallocate(InteractionHash)
201     allocate(InteractionHash(nAtypes,nAtypes))
202 chuckv 561 endif
203 gezelter 571
204     if (.not. allocated(atypeMaxCutoff)) then
205     allocate(atypeMaxCutoff(nAtypes))
206 chuckv 655 else
207     deallocate(atypeMaxCutoff)
208     allocate(atypeMaxCutoff(nAtypes))
209 gezelter 571 endif
210 chuckv 561
211     do i = 1, nAtypes
212     call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
213     call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
214     call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
215     call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
216     call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
217     call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
218     call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
219 chuckv 733 call getElementProperty(atypes, i, "is_SC", i_is_SC)
220 gezelter 117
221 chuckv 561 do j = i, nAtypes
222 chrisfen 532
223 chuckv 561 iHash = 0
224     myRcut = 0.0_dp
225 gezelter 117
226 chuckv 561 call getElementProperty(atypes, j, "is_LennardJones", j_is_LJ)
227     call getElementProperty(atypes, j, "is_Electrostatic", j_is_Elect)
228     call getElementProperty(atypes, j, "is_Sticky", j_is_Sticky)
229     call getElementProperty(atypes, j, "is_StickyPower", j_is_StickyP)
230     call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
231     call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
232     call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
233 chuckv 733 call getElementProperty(atypes, j, "is_SC", j_is_SC)
234 gezelter 117
235 chuckv 561 if (i_is_LJ .and. j_is_LJ) then
236 gezelter 562 iHash = ior(iHash, LJ_PAIR)
237     endif
238    
239     if (i_is_Elect .and. j_is_Elect) then
240     iHash = ior(iHash, ELECTROSTATIC_PAIR)
241     endif
242    
243     if (i_is_Sticky .and. j_is_Sticky) then
244     iHash = ior(iHash, STICKY_PAIR)
245     endif
246 chuckv 561
247 gezelter 562 if (i_is_StickyP .and. j_is_StickyP) then
248     iHash = ior(iHash, STICKYPOWER_PAIR)
249     endif
250 chuckv 561
251 gezelter 562 if (i_is_EAM .and. j_is_EAM) then
252     iHash = ior(iHash, EAM_PAIR)
253 chuckv 561 endif
254    
255 chuckv 733 if (i_is_SC .and. j_is_SC) then
256     iHash = ior(iHash, SC_PAIR)
257     endif
258    
259 chuckv 561 if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
260     if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
261     if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
262 chuckv 1162
263     if ((i_is_EAM.or.i_is_SC).and.(.not.(j_is_EAM.or.j_is_SC))) iHash = ior(iHash, MNM_PAIR)
264     if ((j_is_EAM.or.j_is_SC).and.(.not.(i_is_EAM.or.i_is_SC))) iHash = ior(iHash, MNM_PAIR)
265 chuckv 561
266     if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
267     if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
268     if (i_is_LJ .and. j_is_Shape) iHash = ior(iHash, SHAPE_LJ)
269    
270    
271 chuckv 570 InteractionHash(i,j) = iHash
272     InteractionHash(j,i) = iHash
273 chuckv 561
274     end do
275    
276     end do
277 tim 568
278 gezelter 571 haveInteractionHash = .true.
279     end subroutine createInteractionHash
280 chuckv 561
281 gezelter 762 subroutine createGtypeCutoffMap()
282 gezelter 569
283 gezelter 574 logical :: i_is_LJ
284     logical :: i_is_Elect
285     logical :: i_is_Sticky
286     logical :: i_is_StickyP
287     logical :: i_is_GB
288     logical :: i_is_EAM
289     logical :: i_is_Shape
290 chuckv 831 logical :: i_is_SC
291 gezelter 587 logical :: GtypeFound
292 chuckv 561
293 gezelter 576 integer :: myStatus, nAtypes, i, j, istart, iend, jstart, jend
294 chuckv 652 integer :: n_in_i, me_i, ia, g, atom1, ja, n_in_j,me_j
295 chuckv 589 integer :: nGroupsInRow
296 chuckv 651 integer :: nGroupsInCol
297     integer :: nGroupTypesRow,nGroupTypesCol
298 gezelter 762 real(kind=dp):: thisSigma, bigSigma, thisRcut, tradRcut, tol
299 gezelter 576 real(kind=dp) :: biggestAtypeCutoff
300 gezelter 571
301     if (.not. haveInteractionHash) then
302 gezelter 762 call createInteractionHash()
303 chuckv 567 endif
304 chuckv 589 #ifdef IS_MPI
305     nGroupsInRow = getNgroupsInRow(plan_group_row)
306 chuckv 651 nGroupsInCol = getNgroupsInCol(plan_group_col)
307 chuckv 589 #endif
308 chuckv 563 nAtypes = getSize(atypes)
309 chuckv 599 ! Set all of the initial cutoffs to zero.
310     atypeMaxCutoff = 0.0_dp
311 gezelter 1313 biggestAtypeCutoff = 0.0_dp
312 gezelter 571 do i = 1, nAtypes
313 gezelter 582 if (SimHasAtype(i)) then
314 gezelter 575 call getElementProperty(atypes, i, "is_LennardJones", i_is_LJ)
315     call getElementProperty(atypes, i, "is_Electrostatic", i_is_Elect)
316     call getElementProperty(atypes, i, "is_Sticky", i_is_Sticky)
317     call getElementProperty(atypes, i, "is_StickyPower", i_is_StickyP)
318     call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
319     call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
320     call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
321 chuckv 831 call getElementProperty(atypes, i, "is_SC", i_is_SC)
322 chuckv 599
323 chrisfen 618 if (haveDefaultCutoffs) then
324     atypeMaxCutoff(i) = defaultRcut
325     else
326     if (i_is_LJ) then
327     thisRcut = getSigma(i) * 2.5_dp
328     if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
329     endif
330     if (i_is_Elect) then
331     thisRcut = defaultRcut
332     if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
333     endif
334     if (i_is_Sticky) then
335     thisRcut = getStickyCut(i)
336     if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
337     endif
338     if (i_is_StickyP) then
339     thisRcut = getStickyPowerCut(i)
340     if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
341     endif
342     if (i_is_GB) then
343     thisRcut = getGayBerneCut(i)
344     if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
345     endif
346     if (i_is_EAM) then
347     thisRcut = getEAMCut(i)
348     if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
349     endif
350     if (i_is_Shape) then
351     thisRcut = getShapeCut(i)
352     if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
353     endif
354 chuckv 831 if (i_is_SC) then
355     thisRcut = getSCCut(i)
356     if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
357     endif
358 gezelter 575 endif
359 gezelter 762
360 gezelter 575 if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
361     biggestAtypeCutoff = atypeMaxCutoff(i)
362     endif
363 chrisfen 618
364 gezelter 574 endif
365 gezelter 575 enddo
366 gezelter 581
367 gezelter 575 istart = 1
368 chuckv 651 jstart = 1
369 gezelter 575 #ifdef IS_MPI
370     iend = nGroupsInRow
371 chuckv 651 jend = nGroupsInCol
372 gezelter 575 #else
373     iend = nGroups
374 chuckv 651 jend = nGroups
375 gezelter 575 #endif
376 gezelter 582
377 gezelter 581 !! allocate the groupToGtype and gtypeMaxCutoff here.
378 chuckv 651 if(.not.allocated(groupToGtypeRow)) then
379     ! allocate(groupToGtype(iend))
380     allocate(groupToGtypeRow(iend))
381     else
382     deallocate(groupToGtypeRow)
383     allocate(groupToGtypeRow(iend))
384 chuckv 583 endif
385 chuckv 651 if(.not.allocated(groupMaxCutoffRow)) then
386     allocate(groupMaxCutoffRow(iend))
387     else
388     deallocate(groupMaxCutoffRow)
389     allocate(groupMaxCutoffRow(iend))
390     end if
391    
392     if(.not.allocated(gtypeMaxCutoffRow)) then
393     allocate(gtypeMaxCutoffRow(iend))
394     else
395     deallocate(gtypeMaxCutoffRow)
396     allocate(gtypeMaxCutoffRow(iend))
397     endif
398    
399    
400     #ifdef IS_MPI
401     ! We only allocate new storage if we are in MPI because Ncol /= Nrow
402 chuckv 652 if(.not.associated(groupToGtypeCol)) then
403 chuckv 651 allocate(groupToGtypeCol(jend))
404     else
405     deallocate(groupToGtypeCol)
406     allocate(groupToGtypeCol(jend))
407     end if
408    
409 tim 833 if(.not.associated(groupMaxCutoffCol)) then
410     allocate(groupMaxCutoffCol(jend))
411 chuckv 651 else
412 tim 833 deallocate(groupMaxCutoffCol)
413     allocate(groupMaxCutoffCol(jend))
414 chuckv 651 end if
415 chuckv 652 if(.not.associated(gtypeMaxCutoffCol)) then
416 chuckv 651 allocate(gtypeMaxCutoffCol(jend))
417     else
418     deallocate(gtypeMaxCutoffCol)
419     allocate(gtypeMaxCutoffCol(jend))
420     end if
421    
422     groupMaxCutoffCol = 0.0_dp
423     gtypeMaxCutoffCol = 0.0_dp
424    
425     #endif
426     groupMaxCutoffRow = 0.0_dp
427     gtypeMaxCutoffRow = 0.0_dp
428    
429    
430 gezelter 582 !! first we do a single loop over the cutoff groups to find the
431     !! largest cutoff for any atypes present in this group. We also
432     !! create gtypes at this point.
433    
434 gezelter 960 tol = 1.0e-6_dp
435 chuckv 651 nGroupTypesRow = 0
436 tim 833 nGroupTypesCol = 0
437 gezelter 581 do i = istart, iend
438 gezelter 575 n_in_i = groupStartRow(i+1) - groupStartRow(i)
439 chuckv 651 groupMaxCutoffRow(i) = 0.0_dp
440 gezelter 581 do ia = groupStartRow(i), groupStartRow(i+1)-1
441     atom1 = groupListRow(ia)
442 gezelter 575 #ifdef IS_MPI
443 gezelter 581 me_i = atid_row(atom1)
444 gezelter 575 #else
445 gezelter 581 me_i = atid(atom1)
446     #endif
447 chuckv 651 if (atypeMaxCutoff(me_i).gt.groupMaxCutoffRow(i)) then
448     groupMaxCutoffRow(i)=atypeMaxCutoff(me_i)
449 gezelter 587 endif
450 gezelter 581 enddo
451 chuckv 651 if (nGroupTypesRow.eq.0) then
452     nGroupTypesRow = nGroupTypesRow + 1
453     gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
454     groupToGtypeRow(i) = nGroupTypesRow
455 gezelter 581 else
456 gezelter 587 GtypeFound = .false.
457 chuckv 651 do g = 1, nGroupTypesRow
458     if ( abs(groupMaxCutoffRow(i) - gtypeMaxCutoffRow(g)).lt.tol) then
459     groupToGtypeRow(i) = g
460 gezelter 587 GtypeFound = .true.
461 gezelter 581 endif
462     enddo
463 gezelter 587 if (.not.GtypeFound) then
464 chuckv 651 nGroupTypesRow = nGroupTypesRow + 1
465     gtypeMaxCutoffRow(nGroupTypesRow) = groupMaxCutoffRow(i)
466     groupToGtypeRow(i) = nGroupTypesRow
467 gezelter 587 endif
468 gezelter 581 endif
469 gezelter 587 enddo
470    
471 chuckv 651 #ifdef IS_MPI
472     do j = jstart, jend
473     n_in_j = groupStartCol(j+1) - groupStartCol(j)
474     groupMaxCutoffCol(j) = 0.0_dp
475     do ja = groupStartCol(j), groupStartCol(j+1)-1
476     atom1 = groupListCol(ja)
477    
478     me_j = atid_col(atom1)
479    
480     if (atypeMaxCutoff(me_j).gt.groupMaxCutoffCol(j)) then
481     groupMaxCutoffCol(j)=atypeMaxCutoff(me_j)
482     endif
483     enddo
484    
485     if (nGroupTypesCol.eq.0) then
486     nGroupTypesCol = nGroupTypesCol + 1
487     gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
488     groupToGtypeCol(j) = nGroupTypesCol
489     else
490     GtypeFound = .false.
491     do g = 1, nGroupTypesCol
492     if ( abs(groupMaxCutoffCol(j) - gtypeMaxCutoffCol(g)).lt.tol) then
493     groupToGtypeCol(j) = g
494     GtypeFound = .true.
495     endif
496     enddo
497     if (.not.GtypeFound) then
498     nGroupTypesCol = nGroupTypesCol + 1
499     gtypeMaxCutoffCol(nGroupTypesCol) = groupMaxCutoffCol(j)
500     groupToGtypeCol(j) = nGroupTypesCol
501     endif
502     endif
503     enddo
504    
505     #else
506     ! Set pointers to information we just found
507     nGroupTypesCol = nGroupTypesRow
508     groupToGtypeCol => groupToGtypeRow
509     gtypeMaxCutoffCol => gtypeMaxCutoffRow
510     groupMaxCutoffCol => groupMaxCutoffRow
511     #endif
512    
513 gezelter 581 !! allocate the gtypeCutoffMap here.
514 chuckv 651 allocate(gtypeCutoffMap(nGroupTypesRow,nGroupTypesCol))
515 gezelter 581 !! then we do a double loop over all the group TYPES to find the cutoff
516     !! map between groups of two types
517 chuckv 651 tradRcut = max(maxval(gtypeMaxCutoffRow),maxval(gtypeMaxCutoffCol))
518    
519 gezelter 762 do i = 1, nGroupTypesRow
520 chuckv 651 do j = 1, nGroupTypesCol
521 gezelter 576
522 gezelter 581 select case(cutoffPolicy)
523 gezelter 582 case(TRADITIONAL_CUTOFF_POLICY)
524 chuckv 651 thisRcut = tradRcut
525 gezelter 582 case(MIX_CUTOFF_POLICY)
526 chuckv 651 thisRcut = 0.5_dp * (gtypeMaxCutoffRow(i) + gtypeMaxCutoffCol(j))
527 gezelter 582 case(MAX_CUTOFF_POLICY)
528 chuckv 651 thisRcut = max(gtypeMaxCutoffRow(i), gtypeMaxCutoffCol(j))
529 gezelter 582 case default
530     call handleError("createGtypeCutoffMap", "Unknown Cutoff Policy")
531     return
532     end select
533     gtypeCutoffMap(i,j)%rcut = thisRcut
534 gezelter 762
535     if (thisRcut.gt.largestRcut) largestRcut = thisRcut
536    
537 gezelter 582 gtypeCutoffMap(i,j)%rcutsq = thisRcut*thisRcut
538 gezelter 585
539 gezelter 762 if (.not.haveSkinThickness) then
540     skinThickness = 1.0_dp
541     endif
542    
543     gtypeCutoffMap(i,j)%rlistsq = (thisRcut + skinThickness)**2
544    
545 chrisfen 618 ! sanity check
546    
547     if (haveDefaultCutoffs) then
548     if (abs(gtypeCutoffMap(i,j)%rcut - defaultRcut).gt.0.0001) then
549     call handleError("createGtypeCutoffMap", "user-specified rCut does not match computed group Cutoff")
550     endif
551     endif
552 gezelter 581 enddo
553     enddo
554 gezelter 762
555 chuckv 651 if(allocated(gtypeMaxCutoffRow)) deallocate(gtypeMaxCutoffRow)
556     if(allocated(groupMaxCutoffRow)) deallocate(groupMaxCutoffRow)
557     if(allocated(atypeMaxCutoff)) deallocate(atypeMaxCutoff)
558     #ifdef IS_MPI
559     if(associated(groupMaxCutoffCol)) deallocate(groupMaxCutoffCol)
560     if(associated(gtypeMaxCutoffCol)) deallocate(gtypeMaxCutoffCol)
561     #endif
562     groupMaxCutoffCol => null()
563     gtypeMaxCutoffCol => null()
564    
565 gezelter 581 haveGtypeCutoffMap = .true.
566 chrisfen 596 end subroutine createGtypeCutoffMap
567 chrisfen 578
568 chrisfen 1129 subroutine setCutoffs(defRcut, defRsw, defSP, defSF)
569 chrisfen 596
570 gezelter 762 real(kind=dp),intent(in) :: defRcut, defRsw
571 gezelter 1386 integer, intent(in) :: defSP, defSF
572 gezelter 762 character(len = statusMsgSize) :: errMsg
573     integer :: localError
574    
575 chrisfen 596 defaultRcut = defRcut
576     defaultRsw = defRsw
577 gezelter 1386
578     if (defSP .ne. 0) then
579     defaultDoShiftPot = .true.
580     else
581     defaultDoShiftPot = .false.
582     endif
583     if (defSF .ne. 0) then
584     defaultDoShiftFrc = .true.
585     else
586     defaultDoShiftFrc = .false.
587     endif
588 chrisfen 1129
589 gezelter 762 if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
590 chrisfen 1129 if (defaultDoShiftFrc) then
591     write(errMsg, *) &
592     'cutoffRadius and switchingRadius are set to the', newline &
593 gezelter 1390 // tab, 'same value. OpenMD will use shifted force', newline &
594 chrisfen 1129 // tab, 'potentials instead of switching functions.'
595    
596     call handleInfo("setCutoffs", errMsg)
597     else
598     write(errMsg, *) &
599     'cutoffRadius and switchingRadius are set to the', newline &
600 gezelter 1390 // tab, 'same value. OpenMD will use shifted', newline &
601 chrisfen 1129 // tab, 'potentials instead of switching functions.'
602    
603     call handleInfo("setCutoffs", errMsg)
604    
605     defaultDoShiftPot = .true.
606     endif
607    
608 gezelter 762 endif
609 gezelter 939
610 gezelter 762 localError = 0
611 chrisfen 1129 call setLJDefaultCutoff( defaultRcut, defaultDoShiftPot, &
612     defaultDoShiftFrc )
613 gezelter 813 call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
614 gezelter 938 call setCutoffEAM( defaultRcut )
615     call setCutoffSC( defaultRcut )
616 chuckv 1162 call setMnMDefaultCutoff( defaultRcut, defaultDoShiftPot, &
617     defaultDoShiftFrc )
618 gezelter 939 call set_switch(defaultRsw, defaultRcut)
619 gezelter 889 call setHmatDangerousRcutValue(defaultRcut)
620 gezelter 939
621 chrisfen 618 haveDefaultCutoffs = .true.
622 gezelter 813 haveGtypeCutoffMap = .false.
623 gezelter 939
624 gezelter 762 end subroutine setCutoffs
625 chrisfen 596
626 gezelter 762 subroutine cWasLame()
627    
628     VisitCutoffsAfterComputing = .true.
629     return
630    
631     end subroutine cWasLame
632    
633 chrisfen 596 subroutine setCutoffPolicy(cutPolicy)
634 gezelter 762
635 chrisfen 596 integer, intent(in) :: cutPolicy
636 gezelter 762
637 chrisfen 596 cutoffPolicy = cutPolicy
638 gezelter 762 haveCutoffPolicy = .true.
639 gezelter 813 haveGtypeCutoffMap = .false.
640 gezelter 762
641 gezelter 576 end subroutine setCutoffPolicy
642 gezelter 1126
643 chrisfen 998 subroutine setBoxDipole()
644    
645     do_box_dipole = .true.
646    
647     end subroutine setBoxDipole
648    
649     subroutine getBoxDipole( box_dipole )
650    
651     real(kind=dp), intent(inout), dimension(3) :: box_dipole
652    
653     box_dipole = boxDipole
654    
655     end subroutine getBoxDipole
656    
657 gezelter 762 subroutine setElectrostaticMethod( thisESM )
658    
659     integer, intent(in) :: thisESM
660    
661     electrostaticSummationMethod = thisESM
662     haveElectrostaticSummationMethod = .true.
663 gezelter 574
664 gezelter 762 end subroutine setElectrostaticMethod
665    
666     subroutine setSkinThickness( thisSkin )
667 gezelter 574
668 gezelter 762 real(kind=dp), intent(in) :: thisSkin
669    
670     skinThickness = thisSkin
671 gezelter 813 haveSkinThickness = .true.
672     haveGtypeCutoffMap = .false.
673 gezelter 762
674     end subroutine setSkinThickness
675    
676     subroutine setSimVariables()
677     SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
678     SIM_uses_EAM = SimUsesEAM()
679     SIM_requires_postpair_calc = SimRequiresPostpairCalc()
680     SIM_requires_prepair_calc = SimRequiresPrepairCalc()
681     SIM_uses_PBC = SimUsesPBC()
682 chuckv 841 SIM_uses_SC = SimUsesSC()
683 gezelter 1126 SIM_uses_AtomicVirial = SimUsesAtomicVirial()
684 chrisfen 998
685 gezelter 762 haveSIMvariables = .true.
686    
687     return
688     end subroutine setSimVariables
689 gezelter 117
690     subroutine doReadyCheck(error)
691     integer, intent(out) :: error
692     integer :: myStatus
693    
694     error = 0
695 chrisfen 532
696 gezelter 571 if (.not. haveInteractionHash) then
697 gezelter 762 call createInteractionHash()
698 gezelter 117 endif
699    
700 gezelter 571 if (.not. haveGtypeCutoffMap) then
701 gezelter 762 call createGtypeCutoffMap()
702 gezelter 571 endif
703    
704 gezelter 762 if (VisitCutoffsAfterComputing) then
705 gezelter 939 call set_switch(largestRcut, largestRcut)
706 gezelter 889 call setHmatDangerousRcutValue(largestRcut)
707 gezelter 938 call setCutoffEAM(largestRcut)
708     call setCutoffSC(largestRcut)
709     VisitCutoffsAfterComputing = .false.
710 gezelter 762 endif
711    
712 gezelter 117 if (.not. haveSIMvariables) then
713     call setSimVariables()
714     endif
715    
716     if (.not. haveNeighborList) then
717     write(default_error, *) 'neighbor list has not been initialized in doForces!'
718     error = -1
719     return
720     end if
721 gezelter 939
722 gezelter 117 if (.not. haveSaneForceField) then
723     write(default_error, *) 'Force Field is not sane in doForces!'
724     error = -1
725     return
726     end if
727 gezelter 939
728 gezelter 117 #ifdef IS_MPI
729     if (.not. isMPISimSet()) then
730     write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
731     error = -1
732     return
733     endif
734     #endif
735     return
736     end subroutine doReadyCheck
737    
738 chrisfen 532
739 gezelter 762 subroutine init_FF(thisStat)
740 gezelter 117
741     integer, intent(out) :: thisStat
742     integer :: my_status, nMatches
743     integer, pointer :: MatchList(:) => null()
744    
745     !! assume things are copacetic, unless they aren't
746     thisStat = 0
747    
748     !! init_FF is called *after* all of the atom types have been
749     !! defined in atype_module using the new_atype subroutine.
750     !!
751     !! this will scan through the known atypes and figure out what
752     !! interactions are used by the force field.
753 chrisfen 532
754 gezelter 141 FF_uses_DirectionalAtoms = .false.
755     FF_uses_Dipoles = .false.
756     FF_uses_GayBerne = .false.
757 gezelter 117 FF_uses_EAM = .false.
758 chuckv 834 FF_uses_SC = .false.
759 chrisfen 532
760 gezelter 141 call getMatchingElementList(atypes, "is_Directional", .true., &
761     nMatches, MatchList)
762     if (nMatches .gt. 0) FF_uses_DirectionalAtoms = .true.
763    
764     call getMatchingElementList(atypes, "is_Dipole", .true., &
765     nMatches, MatchList)
766 gezelter 571 if (nMatches .gt. 0) FF_uses_Dipoles = .true.
767 chrisfen 523
768 gezelter 141 call getMatchingElementList(atypes, "is_GayBerne", .true., &
769     nMatches, MatchList)
770 gezelter 571 if (nMatches .gt. 0) FF_uses_GayBerne = .true.
771 chrisfen 532
772 gezelter 117 call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
773     if (nMatches .gt. 0) FF_uses_EAM = .true.
774 chrisfen 532
775 chuckv 834 call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
776     if (nMatches .gt. 0) FF_uses_SC = .true.
777 gezelter 141
778 chuckv 834
779 gezelter 117 haveSaneForceField = .true.
780 chrisfen 532
781 gezelter 117
782     if (.not. haveNeighborList) then
783     !! Create neighbor lists
784     call expandNeighborList(nLocal, my_status)
785     if (my_Status /= 0) then
786     write(default_error,*) "SimSetup: ExpandNeighborList returned error."
787     thisStat = -1
788     return
789     endif
790     haveNeighborList = .true.
791 chrisfen 532 endif
792    
793 gezelter 117 end subroutine init_FF
794    
795 chrisfen 532
796 gezelter 117 !! Does force loop over i,j pairs. Calls do_pair to calculates forces.
797     !------------------------------------------------------------->
798 gezelter 1285 subroutine do_force_loop(q, q_group, A, eFrame, f, t, tau, pot, particle_pot, &
799 gezelter 1464 error)
800 gezelter 117 !! Position array provided by C, dimensioned by getNlocal
801     real ( kind = dp ), dimension(3, nLocal) :: q
802     !! molecular center-of-mass position array
803     real ( kind = dp ), dimension(3, nGroups) :: q_group
804     !! Rotation Matrix for each long range particle in simulation.
805     real( kind = dp), dimension(9, nLocal) :: A
806     !! Unit vectors for dipoles (lab frame)
807 gezelter 246 real( kind = dp ), dimension(9,nLocal) :: eFrame
808 gezelter 117 !! Force array provided by C, dimensioned by getNlocal
809     real ( kind = dp ), dimension(3,nLocal) :: f
810     !! Torsion array provided by C, dimensioned by getNlocal
811     real( kind = dp ), dimension(3,nLocal) :: t
812    
813     !! Stress Tensor
814     real( kind = dp), dimension(9) :: tau
815 gezelter 662 real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
816 chuckv 1245 real( kind = dp ), dimension(nLocal) :: particle_pot
817 gezelter 1464
818 gezelter 117 logical :: in_switching_region
819     #ifdef IS_MPI
820 gezelter 662 real( kind = DP ), dimension(LR_POT_TYPES) :: pot_local
821 gezelter 117 integer :: nAtomsInRow
822     integer :: nAtomsInCol
823     integer :: nprocs
824     integer :: nGroupsInRow
825     integer :: nGroupsInCol
826     #endif
827     integer :: natoms
828     logical :: update_nlist
829     integer :: i, j, jstart, jend, jnab
830     integer :: istart, iend
831     integer :: ia, jb, atom1, atom2
832     integer :: nlist
833 gezelter 1126 real( kind = DP ) :: ratmsq, rgrpsq, rgrp, rag, vpair, vij
834 gezelter 117 real( kind = DP ) :: sw, dswdr, swderiv, mf
835 chrisfen 699 real( kind = DP ) :: rVal
836 gezelter 1126 real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij, fg, dag
837     real(kind=dp) :: rfpot, mu_i
838 gezelter 762 real(kind=dp):: rCut
839 gezelter 1345 integer :: me_i, me_j, n_in_i, n_in_j, iG, j1
840 gezelter 117 logical :: is_dp_i
841     integer :: neighborListSize
842     integer :: listerror, error
843     integer :: localError
844     integer :: propPack_i, propPack_j
845     integer :: loopStart, loopEnd, loop
846 gezelter 1352 integer :: iHash, jHash
847 gezelter 1286 integer :: i1, topoDist
848 chrisfen 532
849 chrisfen 998 !! the variables for the box dipole moment
850     #ifdef IS_MPI
851     integer :: pChgCount_local
852     integer :: nChgCount_local
853     real(kind=dp) :: pChg_local
854     real(kind=dp) :: nChg_local
855     real(kind=dp), dimension(3) :: pChgPos_local
856     real(kind=dp), dimension(3) :: nChgPos_local
857     real(kind=dp), dimension(3) :: dipVec_local
858     #endif
859     integer :: pChgCount
860     integer :: nChgCount
861     real(kind=dp) :: pChg
862     real(kind=dp) :: nChg
863     real(kind=dp) :: chg_value
864     real(kind=dp), dimension(3) :: pChgPos
865     real(kind=dp), dimension(3) :: nChgPos
866     real(kind=dp), dimension(3) :: dipVec
867     real(kind=dp), dimension(3) :: chgVec
868 gezelter 1340 real(kind=dp) :: skch
869 chrisfen 998
870     !! initialize box dipole variables
871     if (do_box_dipole) then
872     #ifdef IS_MPI
873     pChg_local = 0.0_dp
874     nChg_local = 0.0_dp
875     pChgCount_local = 0
876     nChgCount_local = 0
877     do i=1, 3
878     pChgPos_local = 0.0_dp
879     nChgPos_local = 0.0_dp
880     dipVec_local = 0.0_dp
881     enddo
882     #endif
883     pChg = 0.0_dp
884     nChg = 0.0_dp
885     pChgCount = 0
886     nChgCount = 0
887     chg_value = 0.0_dp
888    
889     do i=1, 3
890     pChgPos(i) = 0.0_dp
891     nChgPos(i) = 0.0_dp
892     dipVec(i) = 0.0_dp
893     chgVec(i) = 0.0_dp
894     boxDipole(i) = 0.0_dp
895     enddo
896     endif
897    
898 gezelter 117 !! initialize local variables
899 chrisfen 532
900 gezelter 117 #ifdef IS_MPI
901     pot_local = 0.0_dp
902     nAtomsInRow = getNatomsInRow(plan_atom_row)
903     nAtomsInCol = getNatomsInCol(plan_atom_col)
904     nGroupsInRow = getNgroupsInRow(plan_group_row)
905     nGroupsInCol = getNgroupsInCol(plan_group_col)
906     #else
907     natoms = nlocal
908     #endif
909 chrisfen 532
910 gezelter 117 call doReadyCheck(localError)
911     if ( localError .ne. 0 ) then
912     call handleError("do_force_loop", "Not Initialized")
913     error = -1
914     return
915     end if
916     call zero_work_arrays()
917 chrisfen 532
918 gezelter 117 ! Gather all information needed by all force loops:
919 chrisfen 532
920 gezelter 117 #ifdef IS_MPI
921 chrisfen 532
922 gezelter 117 call gather(q, q_Row, plan_atom_row_3d)
923     call gather(q, q_Col, plan_atom_col_3d)
924    
925     call gather(q_group, q_group_Row, plan_group_row_3d)
926     call gather(q_group, q_group_Col, plan_group_col_3d)
927 chrisfen 532
928 gezelter 141 if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
929 gezelter 246 call gather(eFrame, eFrame_Row, plan_atom_row_rotation)
930     call gather(eFrame, eFrame_Col, plan_atom_col_rotation)
931 chrisfen 532
932 gezelter 117 call gather(A, A_Row, plan_atom_row_rotation)
933     call gather(A, A_Col, plan_atom_col_rotation)
934     endif
935 chrisfen 532
936 gezelter 117 #endif
937 chrisfen 532
938 gezelter 117 !! Begin force loop timing:
939     #ifdef PROFILE
940     call cpu_time(forceTimeInitial)
941     nloops = nloops + 1
942     #endif
943 chrisfen 532
944 gezelter 117 loopEnd = PAIR_LOOP
945     if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
946     loopStart = PREPAIR_LOOP
947     else
948     loopStart = PAIR_LOOP
949     endif
950    
951     do loop = loopStart, loopEnd
952    
953     ! See if we need to update neighbor lists
954     ! (but only on the first time through):
955     if (loop .eq. loopStart) then
956     #ifdef IS_MPI
957 gezelter 762 call checkNeighborList(nGroupsInRow, q_group_row, skinThickness, &
958 chrisfen 532 update_nlist)
959 gezelter 117 #else
960 gezelter 762 call checkNeighborList(nGroups, q_group, skinThickness, &
961 chrisfen 532 update_nlist)
962 gezelter 117 #endif
963     endif
964 chrisfen 532
965 gezelter 117 if (update_nlist) then
966     !! save current configuration and construct neighbor list
967     #ifdef IS_MPI
968     call saveNeighborList(nGroupsInRow, q_group_row)
969     #else
970     call saveNeighborList(nGroups, q_group)
971     #endif
972     neighborListSize = size(list)
973     nlist = 0
974     endif
975 chrisfen 532
976 gezelter 117 istart = 1
977     #ifdef IS_MPI
978     iend = nGroupsInRow
979     #else
980     iend = nGroups - 1
981     #endif
982     outer: do i = istart, iend
983    
984     if (update_nlist) point(i) = nlist + 1
985 chrisfen 532
986 gezelter 117 n_in_i = groupStartRow(i+1) - groupStartRow(i)
987 chrisfen 532
988 gezelter 117 if (update_nlist) then
989     #ifdef IS_MPI
990     jstart = 1
991     jend = nGroupsInCol
992     #else
993     jstart = i+1
994     jend = nGroups
995     #endif
996     else
997     jstart = point(i)
998     jend = point(i+1) - 1
999     ! make sure group i has neighbors
1000     if (jstart .gt. jend) cycle outer
1001     endif
1002 chrisfen 532
1003 gezelter 117 do jnab = jstart, jend
1004     if (update_nlist) then
1005     j = jnab
1006     else
1007     j = list(jnab)
1008     endif
1009    
1010     #ifdef IS_MPI
1011 chuckv 567 me_j = atid_col(j)
1012 gezelter 117 call get_interatomic_vector(q_group_Row(:,i), &
1013     q_group_Col(:,j), d_grp, rgrpsq)
1014     #else
1015 chuckv 567 me_j = atid(j)
1016 gezelter 117 call get_interatomic_vector(q_group(:,i), &
1017     q_group(:,j), d_grp, rgrpsq)
1018 chrisfen 618 #endif
1019 gezelter 117
1020 chuckv 651 if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rListsq) then
1021 gezelter 117 if (update_nlist) then
1022     nlist = nlist + 1
1023 chrisfen 532
1024 gezelter 117 if (nlist > neighborListSize) then
1025     #ifdef IS_MPI
1026     call expandNeighborList(nGroupsInRow, listerror)
1027     #else
1028     call expandNeighborList(nGroups, listerror)
1029     #endif
1030     if (listerror /= 0) then
1031     error = -1
1032     write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
1033     return
1034     end if
1035     neighborListSize = size(list)
1036     endif
1037 chrisfen 532
1038 gezelter 117 list(nlist) = j
1039     endif
1040 gezelter 939
1041 chrisfen 708 if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
1042 chrisfen 532
1043 gezelter 762 rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
1044 chrisfen 708 if (loop .eq. PAIR_LOOP) then
1045 gezelter 960 vij = 0.0_dp
1046 gezelter 938 fij(1) = 0.0_dp
1047     fij(2) = 0.0_dp
1048     fij(3) = 0.0_dp
1049 chrisfen 708 endif
1050    
1051 gezelter 939 call get_switch(rgrpsq, sw, dswdr,rgrp, in_switching_region)
1052 chrisfen 708
1053     n_in_j = groupStartCol(j+1) - groupStartCol(j)
1054    
1055     do ia = groupStartRow(i), groupStartRow(i+1)-1
1056 chrisfen 703
1057 chrisfen 708 atom1 = groupListRow(ia)
1058    
1059     inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
1060    
1061     atom2 = groupListCol(jb)
1062    
1063     if (skipThisPair(atom1, atom2)) cycle inner
1064    
1065     if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1066 gezelter 938 d_atm(1) = d_grp(1)
1067     d_atm(2) = d_grp(2)
1068     d_atm(3) = d_grp(3)
1069 chrisfen 708 ratmsq = rgrpsq
1070     else
1071 gezelter 117 #ifdef IS_MPI
1072 chrisfen 708 call get_interatomic_vector(q_Row(:,atom1), &
1073     q_Col(:,atom2), d_atm, ratmsq)
1074 gezelter 117 #else
1075 chrisfen 708 call get_interatomic_vector(q(:,atom1), &
1076     q(:,atom2), d_atm, ratmsq)
1077 gezelter 117 #endif
1078 gezelter 1286 endif
1079    
1080     topoDist = getTopoDistance(atom1, atom2)
1081    
1082 chrisfen 708 if (loop .eq. PREPAIR_LOOP) then
1083 gezelter 117 #ifdef IS_MPI
1084 chrisfen 708 call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1085 gezelter 1464 rgrpsq, d_grp, rCut, &
1086 chrisfen 708 eFrame, A, f, t, pot_local)
1087 gezelter 117 #else
1088 chrisfen 708 call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
1089 gezelter 1464 rgrpsq, d_grp, rCut, &
1090 chrisfen 708 eFrame, A, f, t, pot)
1091 gezelter 117 #endif
1092 chrisfen 708 else
1093 gezelter 117 #ifdef IS_MPI
1094 chrisfen 708 call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1095 gezelter 1464 eFrame, A, f, t, pot_local, particle_pot, vpair, &
1096 gezelter 1286 fpair, d_grp, rgrp, rCut, topoDist)
1097 chuckv 1245 ! particle_pot will be accumulated from row & column
1098     ! arrays later
1099 gezelter 117 #else
1100 chrisfen 708 call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1101 gezelter 1464 eFrame, A, f, t, pot, particle_pot, vpair, &
1102 gezelter 1286 fpair, d_grp, rgrp, rCut, topoDist)
1103 gezelter 117 #endif
1104 chrisfen 708 vij = vij + vpair
1105 gezelter 938 fij(1) = fij(1) + fpair(1)
1106     fij(2) = fij(2) + fpair(2)
1107     fij(3) = fij(3) + fpair(3)
1108 gezelter 1464 call add_stress_tensor(d_atm, fpair, tau)
1109 chrisfen 708 endif
1110     enddo inner
1111     enddo
1112 gezelter 117
1113 chrisfen 708 if (loop .eq. PAIR_LOOP) then
1114     if (in_switching_region) then
1115     swderiv = vij*dswdr/rgrp
1116 chrisfen 1131 fg = swderiv*d_grp
1117     fij(1) = fij(1) + fg(1)
1118     fij(2) = fij(2) + fg(2)
1119     fij(3) = fij(3) + fg(3)
1120 chrisfen 708
1121 gezelter 1464 if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1122 chrisfen 1131 call add_stress_tensor(d_atm, fg, tau)
1123 gezelter 1464 endif
1124 chrisfen 1131
1125 chrisfen 708 do ia=groupStartRow(i), groupStartRow(i+1)-1
1126     atom1=groupListRow(ia)
1127     mf = mfactRow(atom1)
1128 gezelter 1126 ! fg is the force on atom ia due to cutoff group's
1129     ! presence in switching region
1130     fg = swderiv*d_grp*mf
1131 gezelter 117 #ifdef IS_MPI
1132 gezelter 1126 f_Row(1,atom1) = f_Row(1,atom1) + fg(1)
1133     f_Row(2,atom1) = f_Row(2,atom1) + fg(2)
1134     f_Row(3,atom1) = f_Row(3,atom1) + fg(3)
1135 gezelter 117 #else
1136 gezelter 1126 f(1,atom1) = f(1,atom1) + fg(1)
1137     f(2,atom1) = f(2,atom1) + fg(2)
1138     f(3,atom1) = f(3,atom1) + fg(3)
1139 gezelter 117 #endif
1140 gezelter 1127 if (n_in_i .gt. 1) then
1141 gezelter 1464 if (SIM_uses_AtomicVirial) then
1142     ! find the distance between the atom
1143     ! and the center of the cutoff group:
1144 gezelter 1126 #ifdef IS_MPI
1145 gezelter 1127 call get_interatomic_vector(q_Row(:,atom1), &
1146     q_group_Row(:,i), dag, rag)
1147 gezelter 1126 #else
1148 gezelter 1127 call get_interatomic_vector(q(:,atom1), &
1149     q_group(:,i), dag, rag)
1150 gezelter 1126 #endif
1151 gezelter 1127 call add_stress_tensor(dag,fg,tau)
1152     endif
1153 gezelter 1126 endif
1154 chrisfen 708 enddo
1155    
1156     do jb=groupStartCol(j), groupStartCol(j+1)-1
1157     atom2=groupListCol(jb)
1158     mf = mfactCol(atom2)
1159 gezelter 1126 ! fg is the force on atom jb due to cutoff group's
1160     ! presence in switching region
1161     fg = -swderiv*d_grp*mf
1162 gezelter 117 #ifdef IS_MPI
1163 gezelter 1126 f_Col(1,atom2) = f_Col(1,atom2) + fg(1)
1164     f_Col(2,atom2) = f_Col(2,atom2) + fg(2)
1165     f_Col(3,atom2) = f_Col(3,atom2) + fg(3)
1166 gezelter 117 #else
1167 gezelter 1126 f(1,atom2) = f(1,atom2) + fg(1)
1168     f(2,atom2) = f(2,atom2) + fg(2)
1169     f(3,atom2) = f(3,atom2) + fg(3)
1170 gezelter 117 #endif
1171 gezelter 1127 if (n_in_j .gt. 1) then
1172 gezelter 1464 if (SIM_uses_AtomicVirial) then
1173     ! find the distance between the atom
1174     ! and the center of the cutoff group:
1175 gezelter 1126 #ifdef IS_MPI
1176 gezelter 1127 call get_interatomic_vector(q_Col(:,atom2), &
1177     q_group_Col(:,j), dag, rag)
1178 gezelter 1126 #else
1179 gezelter 1127 call get_interatomic_vector(q(:,atom2), &
1180     q_group(:,j), dag, rag)
1181 gezelter 1126 #endif
1182 gezelter 1464 call add_stress_tensor(dag,fg,tau)
1183 gezelter 1127 endif
1184 gezelter 1464 endif
1185 chrisfen 708 enddo
1186     endif
1187 gezelter 1464 !if (.not.SIM_uses_AtomicVirial) then
1188 gezelter 1174 ! call add_stress_tensor(d_grp, fij, tau)
1189     !endif
1190 gezelter 117 endif
1191     endif
1192 chrisfen 708 endif
1193 gezelter 117 enddo
1194 chrisfen 708
1195 gezelter 117 enddo outer
1196 chrisfen 532
1197 gezelter 117 if (update_nlist) then
1198     #ifdef IS_MPI
1199     point(nGroupsInRow + 1) = nlist + 1
1200     #else
1201     point(nGroups) = nlist + 1
1202     #endif
1203     if (loop .eq. PREPAIR_LOOP) then
1204     ! we just did the neighbor list update on the first
1205     ! pass, so we don't need to do it
1206     ! again on the second pass
1207     update_nlist = .false.
1208     endif
1209     endif
1210 chrisfen 532
1211 gezelter 117 if (loop .eq. PREPAIR_LOOP) then
1212 chuckv 1133 #ifdef IS_MPI
1213 gezelter 1285 call do_preforce(nlocal, pot_local, particle_pot)
1214 chuckv 1133 #else
1215 gezelter 1285 call do_preforce(nlocal, pot, particle_pot)
1216 chuckv 1133 #endif
1217 gezelter 117 endif
1218 chrisfen 532
1219 gezelter 117 enddo
1220 chrisfen 532
1221 gezelter 117 !! Do timing
1222     #ifdef PROFILE
1223     call cpu_time(forceTimeFinal)
1224     forceTime = forceTime + forceTimeFinal - forceTimeInitial
1225     #endif
1226 chrisfen 532
1227 gezelter 117 #ifdef IS_MPI
1228     !!distribute forces
1229 chrisfen 532
1230 gezelter 117 f_temp = 0.0_dp
1231     call scatter(f_Row,f_temp,plan_atom_row_3d)
1232     do i = 1,nlocal
1233     f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
1234     end do
1235 chrisfen 532
1236 gezelter 117 f_temp = 0.0_dp
1237     call scatter(f_Col,f_temp,plan_atom_col_3d)
1238     do i = 1,nlocal
1239     f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
1240     end do
1241 chrisfen 532
1242 gezelter 141 if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
1243 gezelter 117 t_temp = 0.0_dp
1244     call scatter(t_Row,t_temp,plan_atom_row_3d)
1245     do i = 1,nlocal
1246     t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
1247     end do
1248     t_temp = 0.0_dp
1249     call scatter(t_Col,t_temp,plan_atom_col_3d)
1250 chrisfen 532
1251 gezelter 117 do i = 1,nlocal
1252     t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
1253     end do
1254     endif
1255 chrisfen 532
1256 gezelter 1464 ! scatter/gather pot_row into the members of my column
1257     do i = 1,LR_POT_TYPES
1258     call scatter(pot_Row(i,:), pot_Temp(i,:), plan_atom_row)
1259     end do
1260     ! scatter/gather pot_local into all other procs
1261     ! add resultant to get total pot
1262     do i = 1, nlocal
1263     pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES) &
1264     + pot_Temp(1:LR_POT_TYPES,i)
1265     enddo
1266    
1267     do i = 1,LR_POT_TYPES
1268     particle_pot(1:nlocal) = particle_pot(1:nlocal) + pot_Temp(i,1:nlocal)
1269     enddo
1270    
1271     pot_Temp = 0.0_DP
1272    
1273     do i = 1,LR_POT_TYPES
1274     call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1275     end do
1276    
1277     do i = 1, nlocal
1278     pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1279     + pot_Temp(1:LR_POT_TYPES,i)
1280     enddo
1281    
1282     do i = 1,LR_POT_TYPES
1283     particle_pot(1:nlocal) = particle_pot(1:nlocal) + pot_Temp(i,1:nlocal)
1284     enddo
1285    
1286     ppot_Temp = 0.0_DP
1287    
1288     call scatter(ppot_Row(:), ppot_Temp(:), plan_atom_row)
1289     do i = 1, nlocal
1290     particle_pot(i) = particle_pot(i) + ppot_Temp(i)
1291     enddo
1292    
1293     ppot_Temp = 0.0_DP
1294    
1295     call scatter(ppot_Col(:), ppot_Temp(:), plan_atom_col)
1296     do i = 1, nlocal
1297     particle_pot(i) = particle_pot(i) + ppot_Temp(i)
1298     enddo
1299    
1300 gezelter 117 #endif
1301 chrisfen 532
1302 chrisfen 691 if (SIM_requires_postpair_calc) then
1303 chrisfen 695 do i = 1, nlocal
1304    
1305     ! we loop only over the local atoms, so we don't need row and column
1306     ! lookups for the types
1307 gezelter 1346
1308 chrisfen 691 me_i = atid(i)
1309 gezelter 1486
1310 chrisfen 695 ! is the atom electrostatic? See if it would have an
1311     ! electrostatic interaction with itself
1312     iHash = InteractionHash(me_i,me_i)
1313 chrisfen 699
1314 chrisfen 691 if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1315 gezelter 1340
1316     ! loop over the excludes to accumulate charge in the
1317     ! cutoff sphere that we've left out of the normal pair loop
1318     skch = 0.0_dp
1319 gezelter 1345
1320 gezelter 1346 do i1 = 1, nSkipsForLocalAtom(i)
1321     j = skipsForLocalAtom(i, i1)
1322 gezelter 1340 me_j = atid(j)
1323 gezelter 1352 jHash = InteractionHash(me_i,me_j)
1324     if ( iand(jHash, ELECTROSTATIC_PAIR).ne.0 ) then
1325     skch = skch + getCharge(me_j)
1326     endif
1327 gezelter 1340 enddo
1328 gezelter 1346
1329 gezelter 117 #ifdef IS_MPI
1330 gezelter 1464 call self_self(i, eFrame, skch, pot_local(ELECTROSTATIC_POT), t)
1331 gezelter 117 #else
1332 gezelter 1464 call self_self(i, eFrame, skch, pot(ELECTROSTATIC_POT), t)
1333 gezelter 117 #endif
1334 chrisfen 691 endif
1335 chrisfen 699
1336 chrisfen 703
1337 chrisfen 708 if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1338 chrisfen 699
1339 chrisfen 703 ! loop over the excludes to accumulate RF stuff we've
1340     ! left out of the normal pair loop
1341    
1342 gezelter 1346 do i1 = 1, nSkipsForLocalAtom(i)
1343     j = skipsForLocalAtom(i, i1)
1344 chrisfen 703
1345     ! prevent overcounting of the skips
1346     if (i.lt.j) then
1347 gezelter 939 call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq)
1348 gezelter 960 rVal = sqrt(ratmsq)
1349 gezelter 939 call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region)
1350 chrisfen 699 #ifdef IS_MPI
1351 gezelter 1286 call rf_self_excludes(i, j, sw, 1.0_dp, eFrame, d_atm, rVal, &
1352 gezelter 1464 vpair, pot_local(ELECTROSTATIC_POT), f, t)
1353 chrisfen 699 #else
1354 gezelter 1286 call rf_self_excludes(i, j, sw, 1.0_dp, eFrame, d_atm, rVal, &
1355 gezelter 1464 vpair, pot(ELECTROSTATIC_POT), f, t)
1356 chrisfen 699 #endif
1357 chrisfen 703 endif
1358     enddo
1359 chrisfen 708 endif
1360 chrisfen 998
1361     if (do_box_dipole) then
1362     #ifdef IS_MPI
1363     call accumulate_box_dipole(i, eFrame, q(:,i), pChg_local, &
1364     nChg_local, pChgPos_local, nChgPos_local, dipVec_local, &
1365     pChgCount_local, nChgCount_local)
1366     #else
1367     call accumulate_box_dipole(i, eFrame, q(:,i), pChg, nChg, &
1368     pChgPos, nChgPos, dipVec, pChgCount, nChgCount)
1369     #endif
1370     endif
1371 chrisfen 703 enddo
1372 gezelter 117 endif
1373 chrisfen 998
1374 gezelter 117 #ifdef IS_MPI
1375 gezelter 962 #ifdef SINGLE_PRECISION
1376 gezelter 1464 call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_real,mpi_sum, &
1377     mpi_comm_world,mpi_err)
1378 gezelter 962 #else
1379 gezelter 1464 call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision, &
1380     mpi_sum, mpi_comm_world,mpi_err)
1381 gezelter 962 #endif
1382 gezelter 1126
1383 chrisfen 998 if (do_box_dipole) then
1384    
1385     #ifdef SINGLE_PRECISION
1386     call mpi_allreduce(pChg_local, pChg, 1, mpi_real, mpi_sum, &
1387     mpi_comm_world, mpi_err)
1388     call mpi_allreduce(nChg_local, nChg, 1, mpi_real, mpi_sum, &
1389     mpi_comm_world, mpi_err)
1390     call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer, mpi_sum,&
1391     mpi_comm_world, mpi_err)
1392     call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer, mpi_sum,&
1393     mpi_comm_world, mpi_err)
1394     call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_real, mpi_sum, &
1395     mpi_comm_world, mpi_err)
1396     call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_real, mpi_sum, &
1397     mpi_comm_world, mpi_err)
1398     call mpi_allreduce(dipVec_local, dipVec, 3, mpi_real, mpi_sum, &
1399     mpi_comm_world, mpi_err)
1400 gezelter 117 #else
1401 chrisfen 998 call mpi_allreduce(pChg_local, pChg, 1, mpi_double_precision, mpi_sum, &
1402     mpi_comm_world, mpi_err)
1403     call mpi_allreduce(nChg_local, nChg, 1, mpi_double_precision, mpi_sum, &
1404     mpi_comm_world, mpi_err)
1405     call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer,&
1406     mpi_sum, mpi_comm_world, mpi_err)
1407     call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer,&
1408     mpi_sum, mpi_comm_world, mpi_err)
1409     call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_double_precision, &
1410     mpi_sum, mpi_comm_world, mpi_err)
1411     call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_double_precision, &
1412     mpi_sum, mpi_comm_world, mpi_err)
1413     call mpi_allreduce(dipVec_local, dipVec, 3, mpi_double_precision, &
1414     mpi_sum, mpi_comm_world, mpi_err)
1415     #endif
1416    
1417     endif
1418 chrisfen 695
1419 gezelter 117 #endif
1420 chrisfen 998
1421     if (do_box_dipole) then
1422     ! first load the accumulated dipole moment (if dipoles were present)
1423     boxDipole(1) = dipVec(1)
1424     boxDipole(2) = dipVec(2)
1425     boxDipole(3) = dipVec(3)
1426    
1427     ! now include the dipole moment due to charges
1428     ! use the lesser of the positive and negative charge totals
1429     if (nChg .le. pChg) then
1430     chg_value = nChg
1431     else
1432     chg_value = pChg
1433     endif
1434    
1435     ! find the average positions
1436     if (pChgCount .gt. 0 .and. nChgCount .gt. 0) then
1437     pChgPos = pChgPos / pChgCount
1438     nChgPos = nChgPos / nChgCount
1439     endif
1440    
1441     ! dipole is from the negative to the positive (physics notation)
1442     chgVec(1) = pChgPos(1) - nChgPos(1)
1443     chgVec(2) = pChgPos(2) - nChgPos(2)
1444     chgVec(3) = pChgPos(3) - nChgPos(3)
1445    
1446     boxDipole(1) = boxDipole(1) + chgVec(1) * chg_value
1447     boxDipole(2) = boxDipole(2) + chgVec(2) * chg_value
1448     boxDipole(3) = boxDipole(3) + chgVec(3) * chg_value
1449    
1450     endif
1451    
1452 gezelter 117 end subroutine do_force_loop
1453 chrisfen 532
1454 gezelter 1464 subroutine do_pair(i, j, rijsq, d, sw, &
1455 gezelter 1309 eFrame, A, f, t, pot, particle_pot, vpair, &
1456     fpair, d_grp, r_grp, rCut, topoDist)
1457 gezelter 117
1458 chuckv 656 real( kind = dp ) :: vpair, sw
1459 gezelter 662 real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1460 chuckv 1245 real( kind = dp ), dimension(nLocal) :: particle_pot
1461 gezelter 117 real( kind = dp ), dimension(3) :: fpair
1462     real( kind = dp ), dimension(nLocal) :: mfact
1463 gezelter 246 real( kind = dp ), dimension(9,nLocal) :: eFrame
1464 gezelter 117 real( kind = dp ), dimension(9,nLocal) :: A
1465     real( kind = dp ), dimension(3,nLocal) :: f
1466     real( kind = dp ), dimension(3,nLocal) :: t
1467    
1468     integer, intent(in) :: i, j
1469     real ( kind = dp ), intent(inout) :: rijsq
1470 chrisfen 695 real ( kind = dp ), intent(inout) :: r_grp
1471 gezelter 117 real ( kind = dp ), intent(inout) :: d(3)
1472 chrisfen 695 real ( kind = dp ), intent(inout) :: d_grp(3)
1473 gezelter 762 real ( kind = dp ), intent(inout) :: rCut
1474 gezelter 1286 integer, intent(inout) :: topoDist
1475     real ( kind = dp ) :: r, pair_pot, vdwMult, electroMult
1476 gezelter 939 real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1477 gezelter 1386
1478     real( kind = dp), dimension(3) :: f1, t1, t2
1479     real( kind = dp), dimension(9) :: A1, A2, eF1, eF2
1480 chuckv 1388 real( kind = dp) :: dfrhodrho_i, dfrhodrho_j
1481     real( kind = dp) :: rho_i, rho_j
1482     real( kind = dp) :: fshift_i, fshift_j
1483 gezelter 1386 real( kind = dp) :: p_vdw, p_elect, p_hb, p_met
1484     integer :: atid_i, atid_j, id1, id2, idx
1485 gezelter 939 integer :: k
1486 gezelter 117
1487 gezelter 571 integer :: iHash
1488 gezelter 560
1489 chrisfen 942 r = sqrt(rijsq)
1490    
1491 gezelter 960 vpair = 0.0_dp
1492     fpair(1:3) = 0.0_dp
1493 gezelter 117
1494 gezelter 1386 p_vdw = 0.0
1495     p_elect = 0.0
1496     p_hb = 0.0
1497     p_met = 0.0
1498    
1499     f1(1:3) = 0.0
1500    
1501 gezelter 117 #ifdef IS_MPI
1502 gezelter 1386 atid_i = atid_row(i)
1503     atid_j = atid_col(j)
1504 gezelter 117 #else
1505 gezelter 1386 atid_i = atid(i)
1506     atid_j = atid(j)
1507 gezelter 117 #endif
1508 chuckv 1388
1509 gezelter 1386 iHash = InteractionHash(atid_i, atid_j)
1510 cli2 1289
1511 gezelter 1286 vdwMult = vdwScale(topoDist)
1512     electroMult = electrostaticScale(topoDist)
1513 cli2 1289
1514 chrisfen 703 if ( iand(iHash, LJ_PAIR).ne.0 ) then
1515 gezelter 1464 call do_lj_pair(atid_i, atid_j, d, r, rijsq, rcut, sw, vdwMult, vpair, fpair, &
1516     p_vdw, f1)
1517 gezelter 117 endif
1518 chrisfen 532
1519 chrisfen 703 if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1520 gezelter 1510 #ifdef IS_MPI
1521 gezelter 1464 call doElectrostaticPair(atid_i, atid_j, d, r, rijsq, rcut, sw, electroMult, &
1522 gezelter 1510 vpair, fpair, p_elect, eFrame_Row(:,i), eFrame_Col(:,j), &
1523     f1, t_Row(:,i), t_Col(:,j))
1524     #else
1525     call doElectrostaticPair(atid_i, atid_j, d, r, rijsq, rcut, sw, electroMult, &
1526     vpair, fpair, p_elect, eFrame(:,i), eFrame(:,j), f1, t(:,i), t(:,j))
1527     #endif
1528 chrisfen 703 endif
1529    
1530     if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1531 gezelter 1510 #ifdef IS_MPI
1532 gezelter 1464 call do_sticky_pair(atid_i, atid_j, d, r, rijsq, sw, vpair, fpair, &
1533 gezelter 1520 p_hb, A_Row(:,i), A_Col(:,j), f1, t_Row(:,i), t_Col(:,j))
1534 gezelter 1510 #else
1535     call do_sticky_pair(atid_i, atid_j, d, r, rijsq, sw, vpair, fpair, &
1536     p_hb, A(:,i), A(:,j), f1, t(:,i), t(:,j))
1537     #endif
1538 chrisfen 703 endif
1539    
1540     if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1541 gezelter 1510 #ifdef IS_MPI
1542 gezelter 1464 call do_sticky_power_pair(atid_i, atid_j, d, r, rijsq, sw, vpair, fpair, &
1543 gezelter 1520 p_hb, A_Row(:,i), A_Col(:,j), f1, t_Row(:,i), t_Col(:,j))
1544 gezelter 1510 #else
1545     call do_sticky_power_pair(atid_i, atid_j, d, r, rijsq, sw, vpair, fpair, &
1546     p_hb, A(:,i), A(:,j), f1, t(:,i), t(:,j))
1547     #endif
1548 chrisfen 703 endif
1549    
1550     if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1551 gezelter 1510 #ifdef IS_MPI
1552 gezelter 1464 call do_gb_pair(atid_i, atid_j, d, r, rijsq, sw, vdwMult, vpair, fpair, &
1553 gezelter 1520 p_vdw, A_Row(:,i), A_Col(:,j), f1, t_Row(:,i), t_Col(:,j))
1554 gezelter 1510 #else
1555     call do_gb_pair(atid_i, atid_j, d, r, rijsq, sw, vdwMult, vpair, fpair, &
1556     p_vdw, A(:,i), A(:,j), f1, t(:,i), t(:,j))
1557     #endif
1558 chrisfen 703 endif
1559    
1560     if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1561 gezelter 1510 #ifdef IS_MPI
1562 gezelter 1464 call do_gb_pair(atid_i, atid_j, d, r, rijsq, sw, vdwMult, vpair, fpair, &
1563 gezelter 1520 p_vdw, A_Row(:,i), A_Col(:,j), f1, t_Row(:,i), t_Col(:,j))
1564 gezelter 1510 #else
1565     call do_gb_pair(atid_i, atid_j, d, r, rijsq, sw, vdwMult, vpair, fpair, &
1566     p_vdw, A(:,i), A(:,j), f1, t(:,i), t(:,j))
1567     #endif
1568 chrisfen 703 endif
1569 gezelter 1419
1570 chrisfen 703 if ( iand(iHash, SHAPE_PAIR).ne.0 ) then
1571 gezelter 1510 #ifdef IS_MPI
1572 gezelter 1464 call do_shape_pair(atid_i, atid_j, d, r, rijsq, sw, vpair, fpair, &
1573 gezelter 1520 p_vdw, A_Row(:,i), A_Col(:,j), f1, t_Row(:,i), t_Col(:,j))
1574 gezelter 1510 #else
1575     call do_shape_pair(atid_i, atid_j, d, r, rijsq, sw, vpair, fpair, &
1576     p_vdw, A(:,i), A(:,j), f1, t(:,i), t(:,j))
1577     #endif
1578 chrisfen 703 endif
1579    
1580     if ( iand(iHash, SHAPE_LJ).ne.0 ) then
1581 gezelter 1510 #ifdef IS_MPI
1582 gezelter 1464 call do_shape_pair(atid_i, atid_j, d, r, rijsq, sw, vpair, fpair, &
1583 gezelter 1520 p_vdw, A_Row(:,i), A_Col(:,j), f1, t_Row(:,i), t_Col(:,j))
1584 gezelter 1510 #else
1585     call do_shape_pair(atid_i, atid_j, d, r, rijsq, sw, vpair, fpair, &
1586     p_vdw, A(:,i), A(:,j), f1, t(:,i), t(:,j))
1587     #endif
1588 chrisfen 703 endif
1589 chuckv 733
1590 gezelter 1510 if ( iand(iHash, EAM_PAIR).ne.0 ) then
1591     #ifdef IS_MPI
1592 gezelter 1464 call do_eam_pair(atid_i, atid_j, d, r, rijsq, sw, vpair, &
1593 gezelter 1510 fpair, p_met, f1, rho_row(i), rho_col(j), dfrhodrho_row(i), dfrhodrho_col(j), &
1594     fshift_i, fshift_j)
1595     #else
1596     call do_eam_pair(atid_i, atid_j, d, r, rijsq, sw, vpair, &
1597     fpair, p_met, f1, rho(i), rho(j), dfrhodrho(i), dfrhodrho(j), fshift_i, fshift_j)
1598     #endif
1599 gezelter 1419 endif
1600    
1601 chuckv 733 if ( iand(iHash, SC_PAIR).ne.0 ) then
1602 gezelter 1510 #ifdef IS_MPI
1603 gezelter 1464 call do_SC_pair(atid_i, atid_j, d, r, rijsq, sw, vpair, &
1604 gezelter 1510 fpair, p_met, f1, rho_row(i), rho_col(j), dfrhodrho_row(i), dfrhodrho_col(j), &
1605     fshift_i, fshift_j)
1606     #else
1607     call do_SC_pair(atid_i, atid_j, d, r, rijsq, sw, vpair, &
1608     fpair, p_met, f1, rho(i), rho(j), dfrhodrho(i), dfrhodrho(j), fshift_i, fshift_j)
1609     #endif
1610 chuckv 733 endif
1611 chrisfen 703
1612 gezelter 1174 if ( iand(iHash, MNM_PAIR).ne.0 ) then
1613 gezelter 1510 #ifdef IS_MPI
1614 gezelter 1464 call do_mnm_pair(atid_i, atid_j, d, r, rijsq, rcut, sw, vdwMult, vpair, fpair, &
1615 gezelter 1510 p_vdw, A_Row(:,i), A_Col(:,j), f1, t_Row(:,i), t_Col(:,j))
1616     #else
1617     call do_mnm_pair(atid_i, atid_j, d, r, rijsq, rcut, sw, vdwMult, vpair, fpair, &
1618     p_vdw, A(:,i), A(:,j), f1, t(:,i), t(:,j))
1619     #endif
1620 gezelter 1174 endif
1621 gezelter 1386
1622    
1623     #ifdef IS_MPI
1624     id1 = AtomRowToGlobal(i)
1625     id2 = AtomColToGlobal(j)
1626    
1627     pot_row(VDW_POT,i) = pot_row(VDW_POT,i) + 0.5*p_vdw
1628     pot_col(VDW_POT,j) = pot_col(VDW_POT,j) + 0.5*p_vdw
1629     pot_row(ELECTROSTATIC_POT,i) = pot_row(ELECTROSTATIC_POT,i) + 0.5*p_elect
1630     pot_col(ELECTROSTATIC_POT,j) = pot_col(ELECTROSTATIC_POT,j) + 0.5*p_elect
1631     pot_row(HB_POT,i) = pot_row(HB_POT,i) + 0.5*p_hb
1632     pot_col(HB_POT,j) = pot_col(HB_POT,j) + 0.5*p_hb
1633     pot_Row(METALLIC_POT,i) = pot_Row(METALLIC_POT,i) + 0.5*p_met
1634     pot_Col(METALLIC_POT,j) = pot_Col(METALLIC_POT,j) + 0.5*p_met
1635    
1636     do idx = 1, 3
1637     f_Row(idx,i) = f_Row(idx,i) + f1(idx)
1638     f_Col(idx,j) = f_Col(idx,j) - f1(idx)
1639     enddo
1640 chuckv 1388 ! particle_pot is the difference between the full potential
1641     ! and the full potential without the presence of a particular
1642     ! particle (atom1).
1643     !
1644     ! This reduces the density at other particle locations, so
1645     ! we need to recompute the density at atom2 assuming atom1
1646     ! didn't contribute. This then requires recomputing the
1647     ! density functional for atom2 as well.
1648     !
1649     ! Most of the particle_pot heavy lifting comes from the
1650     ! pair interaction, and will be handled by vpair. Parallel version.
1651    
1652 gezelter 1390 if ( (iand(iHash, EAM_PAIR).ne.0) .or. (iand(iHash, SC_PAIR).ne.0) ) then
1653 chuckv 1388 ppot_row(i) = ppot_row(i) - frho_row(j) + fshift_j
1654     ppot_col(j) = ppot_col(j) - frho_col(i) + fshift_i
1655     end if
1656    
1657 gezelter 1386 #else
1658     id1 = i
1659     id2 = j
1660    
1661     pot(VDW_POT) = pot(VDW_POT) + p_vdw
1662     pot(ELECTROSTATIC_POT) = pot(ELECTROSTATIC_POT) + p_elect
1663     pot(HB_POT) = pot(HB_POT) + p_hb
1664     pot(METALLIC_POT) = pot(METALLIC_POT) + p_met
1665    
1666     do idx = 1, 3
1667     f(idx,i) = f(idx,i) + f1(idx)
1668     f(idx,j) = f(idx,j) - f1(idx)
1669     enddo
1670 chuckv 1388 ! particle_pot is the difference between the full potential
1671     ! and the full potential without the presence of a particular
1672     ! particle (atom1).
1673     !
1674     ! This reduces the density at other particle locations, so
1675     ! we need to recompute the density at atom2 assuming atom1
1676     ! didn't contribute. This then requires recomputing the
1677     ! density functional for atom2 as well.
1678     !
1679     ! Most of the particle_pot heavy lifting comes from the
1680     ! pair interaction, and will be handled by vpair. NonParallel version.
1681 gezelter 1390
1682     if ( (iand(iHash, EAM_PAIR).ne.0) .or. (iand(iHash, SC_PAIR).ne.0) ) then
1683 chuckv 1388 particle_pot(i) = particle_pot(i) - frho(j) + fshift_j
1684     particle_pot(j) = particle_pot(j) - frho(i) + fshift_i
1685     end if
1686    
1687    
1688 gezelter 1386 #endif
1689    
1690     if (molMembershipList(id1) .ne. molMembershipList(id2)) then
1691    
1692     fpair(1) = fpair(1) + f1(1)
1693     fpair(2) = fpair(2) + f1(2)
1694     fpair(3) = fpair(3) + f1(3)
1695    
1696     endif
1697    
1698    
1699 gezelter 1313 !!$
1700     !!$ particle_pot(i) = particle_pot(i) + vpair*sw
1701     !!$ particle_pot(j) = particle_pot(j) + vpair*sw
1702 gezelter 1174
1703 gezelter 117 end subroutine do_pair
1704    
1705 gezelter 762 subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
1706 gezelter 1464 eFrame, A, f, t, pot)
1707 gezelter 1390
1708 chuckv 656 real( kind = dp ) :: sw
1709 gezelter 662 real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1710 chrisfen 532 real( kind = dp ), dimension(9,nLocal) :: eFrame
1711     real (kind=dp), dimension(9,nLocal) :: A
1712     real (kind=dp), dimension(3,nLocal) :: f
1713     real (kind=dp), dimension(3,nLocal) :: t
1714 gezelter 1390
1715 chrisfen 532 integer, intent(in) :: i, j
1716 gezelter 762 real ( kind = dp ), intent(inout) :: rijsq, rcijsq, rCut
1717 chrisfen 532 real ( kind = dp ) :: r, rc
1718     real ( kind = dp ), intent(inout) :: d(3), dc(3)
1719 chuckv 1389 real ( kind = dp ) :: rho_i_at_j, rho_j_at_i
1720 gezelter 1386 integer :: atid_i, atid_j, iHash
1721 gezelter 1390
1722 chrisfen 942 r = sqrt(rijsq)
1723    
1724 gezelter 117 #ifdef IS_MPI
1725 gezelter 1386 atid_i = atid_row(i)
1726     atid_j = atid_col(j)
1727 gezelter 117 #else
1728 gezelter 1386 atid_i = atid(i)
1729     atid_j = atid(j)
1730 gezelter 117 #endif
1731 chuckv 1388 rho_i_at_j = 0.0_dp
1732     rho_j_at_i = 0.0_dp
1733    
1734 gezelter 1386 iHash = InteractionHash(atid_i, atid_j)
1735 chrisfen 532
1736 gezelter 571 if ( iand(iHash, EAM_PAIR).ne.0 ) then
1737 gezelter 1464 call calc_EAM_prepair_rho(atid_i, atid_j, d, r, rijsq, rho_i_at_j, rho_j_at_i)
1738 chrisfen 532 endif
1739 gezelter 1390
1740 chuckv 733 if ( iand(iHash, SC_PAIR).ne.0 ) then
1741 gezelter 1464 call calc_SC_prepair_rho(atid_i, atid_j, d, r, rijsq, rho_i_at_j, rho_j_at_i)
1742 chuckv 733 endif
1743 chuckv 1388
1744 gezelter 1390 if ( iand(iHash, EAM_PAIR).ne.0 .or. iand(iHash, SC_PAIR).ne.0 ) then
1745 chuckv 1388 #ifdef IS_MPI
1746     rho_col(j) = rho_col(j) + rho_i_at_j
1747     rho_row(i) = rho_row(i) + rho_j_at_i
1748     #else
1749     rho(j) = rho(j) + rho_i_at_j
1750     rho(i) = rho(i) + rho_j_at_i
1751     #endif
1752     endif
1753 gezelter 560
1754 chrisfen 532 end subroutine do_prepair
1755    
1756    
1757 gezelter 1285 subroutine do_preforce(nlocal, pot, particle_pot)
1758 chrisfen 532 integer :: nlocal
1759 gezelter 662 real( kind = dp ),dimension(LR_POT_TYPES) :: pot
1760 gezelter 1285 real( kind = dp ),dimension(nlocal) :: particle_pot
1761 chuckv 1388 integer :: sc_err = 0
1762 chrisfen 532
1763 chuckv 1388 #ifdef IS_MPI
1764 chuckv 1389 if ((FF_uses_EAM .and. SIM_uses_EAM) .or. (FF_uses_SC .and. SIM_uses_SC)) then
1765 chuckv 1388 call scatter(rho_row,rho,plan_atom_row,sc_err)
1766     if (sc_err /= 0 ) then
1767     call handleError("do_preforce()", "Error scattering rho_row into rho")
1768     endif
1769     call scatter(rho_col,rho_tmp,plan_atom_col,sc_err)
1770     if (sc_err /= 0 ) then
1771     call handleError("do_preforce()", "Error scattering rho_col into rho")
1772     endif
1773     rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal)
1774     end if
1775     #endif
1776    
1777    
1778    
1779 chrisfen 532 if (FF_uses_EAM .and. SIM_uses_EAM) then
1780 gezelter 1285 call calc_EAM_preforce_Frho(nlocal, pot(METALLIC_POT), particle_pot)
1781 chrisfen 532 endif
1782 chuckv 733 if (FF_uses_SC .and. SIM_uses_SC) then
1783 gezelter 1285 call calc_SC_preforce_Frho(nlocal, pot(METALLIC_POT), particle_pot)
1784 chuckv 733 endif
1785 chuckv 1388
1786     #ifdef IS_MPI
1787 chuckv 1389 if ((FF_uses_EAM .and. SIM_uses_EAM) .or. (FF_uses_SC .and. SIM_uses_SC)) then
1788 chuckv 1388 !! communicate f(rho) and derivatives back into row and column arrays
1789     call gather(frho,frho_row,plan_atom_row, sc_err)
1790     if (sc_err /= 0) then
1791     call handleError("do_preforce()","MPI gather frho_row failure")
1792     endif
1793     call gather(dfrhodrho,dfrhodrho_row,plan_atom_row, sc_err)
1794     if (sc_err /= 0) then
1795     call handleError("do_preforce()","MPI gather dfrhodrho_row failure")
1796     endif
1797     call gather(frho,frho_col,plan_atom_col, sc_err)
1798     if (sc_err /= 0) then
1799     call handleError("do_preforce()","MPI gather frho_col failure")
1800     endif
1801     call gather(dfrhodrho,dfrhodrho_col,plan_atom_col, sc_err)
1802     if (sc_err /= 0) then
1803     call handleError("do_preforce()","MPI gather dfrhodrho_col failure")
1804     endif
1805     end if
1806     #endif
1807    
1808 chrisfen 532 end subroutine do_preforce
1809    
1810    
1811     subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1812    
1813     real (kind = dp), dimension(3) :: q_i
1814     real (kind = dp), dimension(3) :: q_j
1815     real ( kind = dp ), intent(out) :: r_sq
1816     real( kind = dp ) :: d(3), scaled(3)
1817 chuckv 1507 real(kind=dp)::t
1818 chrisfen 532 integer i
1819    
1820 gezelter 938 d(1) = q_j(1) - q_i(1)
1821     d(2) = q_j(2) - q_i(2)
1822     d(3) = q_j(3) - q_i(3)
1823 chrisfen 532
1824     ! Wrap back into periodic box if necessary
1825     if ( SIM_uses_PBC ) then
1826    
1827     if( .not.boxIsOrthorhombic ) then
1828     ! calc the scaled coordinates.
1829 gezelter 1508 ! unwrap the matmul and do things explicitly
1830 gezelter 939 ! scaled = matmul(HmatInv, d)
1831 chrisfen 532
1832 gezelter 939 scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1833     scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1834     scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1835    
1836 gezelter 1508 ! wrap the scaled coordinates (but don't use anint for speed)
1837 chrisfen 532
1838 chuckv 1507 t = scaled(1)
1839 gezelter 1517 if (t .gt. 0.0) then
1840 chuckv 1507 scaled(1) = t - floor(t + 0.5)
1841     else
1842 gezelter 1516 scaled(1) = t - ceiling(t - 0.5)
1843 chuckv 1507 endif
1844 chrisfen 532
1845 chuckv 1507 t = scaled(2)
1846 gezelter 1517 if (t .gt. 0.0) then
1847 chuckv 1507 scaled(2) = t - floor(t + 0.5)
1848     else
1849 gezelter 1516 scaled(2) = t - ceiling(t - 0.5)
1850 chuckv 1507 endif
1851    
1852     t = scaled(3)
1853 gezelter 1517 if (t .gt. 0.0) then
1854 chuckv 1507 scaled(3) = t - floor(t + 0.5)
1855     else
1856 gezelter 1516 scaled(3) = t - ceiling(t - 0.5)
1857 chuckv 1507 endif
1858    
1859 chrisfen 532 ! calc the wrapped real coordinates from the wrapped scaled
1860     ! coordinates
1861 gezelter 939 ! d = matmul(Hmat,scaled)
1862     d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1863     d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1864     d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1865 chrisfen 532
1866     else
1867 chuckv 1507 ! calc the scaled coordinates
1868 gezelter 938 scaled(1) = d(1) * HmatInv(1,1)
1869     scaled(2) = d(2) * HmatInv(2,2)
1870     scaled(3) = d(3) * HmatInv(3,3)
1871    
1872     ! wrap the scaled coordinates
1873    
1874 chuckv 1507 t = scaled(1)
1875 gezelter 1517 if (t .gt. 0.0) then
1876 chuckv 1507 scaled(1) = t - floor(t + 0.5)
1877     else
1878 gezelter 1516 scaled(1) = t - ceiling(t - 0.5)
1879 chuckv 1507 endif
1880 chrisfen 532
1881 chuckv 1507 t = scaled(2)
1882 gezelter 1517 if (t .gt. 0.0) then
1883 chuckv 1507 scaled(2) = t - floor(t + 0.5)
1884     else
1885 gezelter 1516 scaled(2) = t - ceiling(t - 0.5)
1886 chuckv 1507 endif
1887    
1888     t = scaled(3)
1889 gezelter 1517 if (t .gt. 0.0) then
1890 chuckv 1507 scaled(3) = t - floor(t + 0.5)
1891     else
1892 gezelter 1516 scaled(3) = t - ceiling(t - 0.5)
1893 chuckv 1507 endif
1894    
1895 gezelter 938 ! calc the wrapped real coordinates from the wrapped scaled
1896     ! coordinates
1897 chrisfen 532
1898 gezelter 938 d(1) = scaled(1)*Hmat(1,1)
1899     d(2) = scaled(2)*Hmat(2,2)
1900     d(3) = scaled(3)*Hmat(3,3)
1901 chrisfen 532
1902     endif
1903    
1904     endif
1905    
1906 gezelter 938 r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1907 chrisfen 532
1908     end subroutine get_interatomic_vector
1909    
1910     subroutine zero_work_arrays()
1911    
1912 gezelter 117 #ifdef IS_MPI
1913    
1914 chrisfen 532 q_Row = 0.0_dp
1915     q_Col = 0.0_dp
1916    
1917     q_group_Row = 0.0_dp
1918     q_group_Col = 0.0_dp
1919    
1920     eFrame_Row = 0.0_dp
1921     eFrame_Col = 0.0_dp
1922    
1923     A_Row = 0.0_dp
1924     A_Col = 0.0_dp
1925    
1926     f_Row = 0.0_dp
1927     f_Col = 0.0_dp
1928     f_Temp = 0.0_dp
1929    
1930     t_Row = 0.0_dp
1931     t_Col = 0.0_dp
1932     t_Temp = 0.0_dp
1933    
1934     pot_Row = 0.0_dp
1935     pot_Col = 0.0_dp
1936     pot_Temp = 0.0_dp
1937 gezelter 1309 ppot_Temp = 0.0_dp
1938 chrisfen 532
1939 chuckv 1388 frho_row = 0.0_dp
1940     frho_col = 0.0_dp
1941     rho_row = 0.0_dp
1942     rho_col = 0.0_dp
1943     rho_tmp = 0.0_dp
1944     dfrhodrho_row = 0.0_dp
1945     dfrhodrho_col = 0.0_dp
1946    
1947 gezelter 117 #endif
1948 chuckv 1388 rho = 0.0_dp
1949     frho = 0.0_dp
1950     dfrhodrho = 0.0_dp
1951 chrisfen 532
1952     end subroutine zero_work_arrays
1953    
1954     function skipThisPair(atom1, atom2) result(skip_it)
1955     integer, intent(in) :: atom1
1956     integer, intent(in), optional :: atom2
1957     logical :: skip_it
1958     integer :: unique_id_1, unique_id_2
1959     integer :: me_i,me_j
1960     integer :: i
1961    
1962     skip_it = .false.
1963    
1964     !! there are a number of reasons to skip a pair or a particle
1965     !! mostly we do this to exclude atoms who are involved in short
1966     !! range interactions (bonds, bends, torsions), but we also need
1967     !! to exclude some overcounted interactions that result from
1968     !! the parallel decomposition
1969    
1970 gezelter 117 #ifdef IS_MPI
1971 chrisfen 532 !! in MPI, we have to look up the unique IDs for each atom
1972     unique_id_1 = AtomRowToGlobal(atom1)
1973     unique_id_2 = AtomColToGlobal(atom2)
1974     !! this situation should only arise in MPI simulations
1975     if (unique_id_1 == unique_id_2) then
1976     skip_it = .true.
1977     return
1978     end if
1979    
1980     !! this prevents us from doing the pair on multiple processors
1981     if (unique_id_1 < unique_id_2) then
1982     if (mod(unique_id_1 + unique_id_2,2) == 0) then
1983     skip_it = .true.
1984     return
1985     endif
1986     else
1987     if (mod(unique_id_1 + unique_id_2,2) == 1) then
1988     skip_it = .true.
1989     return
1990     endif
1991     endif
1992 gezelter 1286 #else
1993     !! in the normal loop, the atom numbers are unique
1994     unique_id_1 = atom1
1995     unique_id_2 = atom2
1996 gezelter 117 #endif
1997 gezelter 1346
1998     #ifdef IS_MPI
1999     do i = 1, nSkipsForRowAtom(atom1)
2000     if (skipsForRowAtom(atom1, i) .eq. unique_id_2) then
2001 chrisfen 532 skip_it = .true.
2002     return
2003     endif
2004     end do
2005 gezelter 1346 #else
2006     do i = 1, nSkipsForLocalAtom(atom1)
2007     if (skipsForLocalAtom(atom1, i) .eq. unique_id_2) then
2008     skip_it = .true.
2009     return
2010     endif
2011     end do
2012     #endif
2013 chrisfen 532
2014     return
2015     end function skipThisPair
2016    
2017 gezelter 1286 function getTopoDistance(atom1, atom2) result(topoDist)
2018     integer, intent(in) :: atom1
2019     integer, intent(in) :: atom2
2020     integer :: topoDist
2021     integer :: unique_id_2
2022     integer :: i
2023    
2024     #ifdef IS_MPI
2025     unique_id_2 = AtomColToGlobal(atom2)
2026     #else
2027     unique_id_2 = atom2
2028     #endif
2029    
2030     ! zero is default for unconnected (i.e. normal) pair interactions
2031    
2032     topoDist = 0
2033    
2034     do i = 1, nTopoPairsForAtom(atom1)
2035     if (toposForAtom(atom1, i) .eq. unique_id_2) then
2036     topoDist = topoDistance(atom1, i)
2037     return
2038     endif
2039     end do
2040    
2041     return
2042     end function getTopoDistance
2043    
2044 chrisfen 532 function FF_UsesDirectionalAtoms() result(doesit)
2045     logical :: doesit
2046 gezelter 571 doesit = FF_uses_DirectionalAtoms
2047 chrisfen 532 end function FF_UsesDirectionalAtoms
2048    
2049     function FF_RequiresPrepairCalc() result(doesit)
2050     logical :: doesit
2051 chuckv 1162 doesit = FF_uses_EAM .or. FF_uses_SC
2052 chrisfen 532 end function FF_RequiresPrepairCalc
2053    
2054 gezelter 117 #ifdef PROFILE
2055 chrisfen 532 function getforcetime() result(totalforcetime)
2056     real(kind=dp) :: totalforcetime
2057     totalforcetime = forcetime
2058     end function getforcetime
2059 gezelter 117 #endif
2060    
2061 chrisfen 532 !! This cleans componets of force arrays belonging only to fortran
2062    
2063 gezelter 1126 subroutine add_stress_tensor(dpair, fpair, tau)
2064 chrisfen 532
2065     real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
2066 gezelter 1126 real( kind = dp ), dimension(9), intent(inout) :: tau
2067 chrisfen 532
2068     ! because the d vector is the rj - ri vector, and
2069     ! because fx, fy, fz are the force on atom i, we need a
2070     ! negative sign here:
2071    
2072 gezelter 1126 tau(1) = tau(1) - dpair(1) * fpair(1)
2073     tau(2) = tau(2) - dpair(1) * fpair(2)
2074     tau(3) = tau(3) - dpair(1) * fpair(3)
2075     tau(4) = tau(4) - dpair(2) * fpair(1)
2076     tau(5) = tau(5) - dpair(2) * fpair(2)
2077     tau(6) = tau(6) - dpair(2) * fpair(3)
2078     tau(7) = tau(7) - dpair(3) * fpair(1)
2079     tau(8) = tau(8) - dpair(3) * fpair(2)
2080     tau(9) = tau(9) - dpair(3) * fpair(3)
2081 chrisfen 532
2082     end subroutine add_stress_tensor
2083    
2084 gezelter 117 end module doForces

Properties

Name Value
svn:keywords Author Id Revision Date