ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/UseTheForce/doForces.F90
Revision: 1503
Committed: Sat Oct 2 19:54:41 2010 UTC (14 years, 7 months ago) by gezelter
File size: 55182 byte(s)
Log Message:
Changes to remove more of the low level stuff from the fortran side.

File Contents

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

Properties

Name Value
svn:keywords Author Id Revision Date