ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/do_Forces.F90
(Generate patch)

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 1131 by tim, Thu Apr 22 21:33:55 2004 UTC vs.
Revision 1199 by gezelter, Thu May 27 15:21:20 2004 UTC

# Line 4 | Line 4
4  
5   !! @author Charles F. Vardeman II
6   !! @author Matthew Meineke
7 < !! @version $Id: do_Forces.F90,v 1.50 2004-04-22 21:33:55 tim Exp $, $Date: 2004-04-22 21:33:55 $, $Name: not supported by cvs2svn $, $Revision: 1.50 $
7 > !! @version $Id: do_Forces.F90,v 1.64 2004-05-27 15:21:20 gezelter Exp $, $Date: 2004-05-27 15:21:20 $, $Name: not supported by cvs2svn $, $Revision: 1.64 $
8  
9   module do_Forces
10    use force_globals
11    use simulation
12    use definitions
13    use atype_module
14 +  use switcheroo
15    use neighborLists  
16    use lj
17    use sticky_pair
# Line 30 | Line 31 | module do_Forces
31  
32   #define __FORTRAN90
33   #include "fForceField.h"
34 + #include "fSwitchingFunction.h"
35 +
36 +  INTEGER, PARAMETER:: PREPAIR_LOOP = 1
37 +  INTEGER, PARAMETER:: PAIR_LOOP    = 2
38  
39    logical, save :: haveRlist = .false.
40    logical, save :: haveNeighborList = .false.
# Line 55 | Line 60 | module do_Forces
60    logical, save :: SIM_requires_prepair_calc
61    logical, save :: SIM_uses_directional_atoms
62    logical, save :: SIM_uses_PBC
63 +  logical, save :: SIM_uses_molecular_cutoffs
64  
65    real(kind=dp), save :: rlist, rlistsq
66  
67    public :: init_FF
68    public :: do_force_loop
69    public :: setRlistDF
64 –
70  
71   #ifdef PROFILE
72    public :: getforcetime
# Line 160 | Line 165 | contains
165      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
166      SIM_uses_directional_atoms = SimUsesDirectionalAtoms()
167      SIM_uses_PBC = SimUsesPBC()
168 +    !SIM_uses_molecular_cutoffs = SimUsesMolecularCutoffs()
169  
170      haveSIMvariables = .true.
171  
# Line 356 | Line 362 | contains
362         endif
363         haveNeighborList = .true.
364      endif
365 +
366      
367 +    
368    end subroutine init_FF
369    
370  
371    !! Does force loop over i,j pairs. Calls do_pair to calculates forces.
372    !------------------------------------------------------------->
373 <  subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
374 <       error)
373 >  subroutine do_force_loop(q, q_group, A, u_l, f, t, tau, pot, &
374 >       do_pot_c, do_stress_c, error)
375      !! Position array provided by C, dimensioned by getNlocal
376 <    real ( kind = dp ), dimension(3,nLocal) :: q
376 >    real ( kind = dp ), dimension(3, nLocal) :: q
377 >    !! molecular center-of-mass position array
378 >    real ( kind = dp ), dimension(3, nGroups) :: q_group
379      !! Rotation Matrix for each long range particle in simulation.
380 <    real( kind = dp), dimension(9,nLocal) :: A    
380 >    real( kind = dp), dimension(9, nLocal) :: A    
381      !! Unit vectors for dipoles (lab frame)
382      real( kind = dp ), dimension(3,nLocal) :: u_l
383      !! Force array provided by C, dimensioned by getNlocal
# Line 381 | Line 391 | contains
391      logical ( kind = 2) :: do_pot_c, do_stress_c
392      logical :: do_pot
393      logical :: do_stress
394 +    logical :: in_switching_region
395   #ifdef IS_MPI
396      real( kind = DP ) :: pot_local
397 <    integer :: nrow
398 <    integer :: ncol
397 >    integer :: nAtomsInRow
398 >    integer :: nAtomsInCol
399      integer :: nprocs
400 +    integer :: nGroupsInRow
401 +    integer :: nGroupsInCol
402   #endif
403      integer :: natoms    
404      logical :: update_nlist  
405 <    integer :: i, j, jbeg, jend, jnab
405 >    integer :: i, j, jstart, jend, jnab
406 >    integer :: istart, iend
407 >    integer :: ia, jb, atom1, atom2
408      integer :: nlist
409 <    real( kind = DP ) ::  rijsq
410 <    real(kind=dp),dimension(3) :: d
409 >    real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
410 >    real( kind = DP ) :: sw, dswdr, swderiv, mf
411 >    real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
412      real(kind=dp) :: rfpot, mu_i, virial
413 <    integer :: me_i, me_j
413 >    integer :: me_i, me_j, n_in_i, n_in_j
414      logical :: is_dp_i
415      integer :: neighborListSize
416      integer :: listerror, error
417      integer :: localError
418      integer :: propPack_i, propPack_j
419 <
404 <    real(kind=dp) :: listSkin = 1.0  
419 >    integer :: loopStart, loopEnd, loop
420  
421 +    real(kind=dp) :: listSkin = 1.0  
422 +    
423      !! initialize local variables  
424 <
424 >    
425   #ifdef IS_MPI
426      pot_local = 0.0_dp
427 <    nrow   = getNrow(plan_row)
428 <    ncol   = getNcol(plan_col)
427 >    nAtomsInRow   = getNatomsInRow(plan_atom_row)
428 >    nAtomsInCol   = getNatomsInCol(plan_atom_col)
429 >    nGroupsInRow  = getNgroupsInRow(plan_group_row)
430 >    nGroupsInCol  = getNgroupsInCol(plan_group_col)
431   #else
432      natoms = nlocal
433   #endif
434 <
434 >    
435      call doReadyCheck(localError)
436      if ( localError .ne. 0 ) then
437         call handleError("do_force_loop", "Not Initialized")
# Line 420 | Line 439 | contains
439         return
440      end if
441      call zero_work_arrays()
442 <
424 <
442 >        
443      do_pot = do_pot_c
444      do_stress = do_stress_c
445 <
445 >    
446      ! Gather all information needed by all force loops:
447      
448   #ifdef IS_MPI    
449 +    
450 +    call gather(q, q_Row, plan_atom_row_3d)
451 +    call gather(q, q_Col, plan_atom_col_3d)
452  
453 <    call gather(q,q_Row,plan_row3d)
454 <    call gather(q,q_Col,plan_col3d)
453 >    call gather(q_group, q_group_Row, plan_group_row_3d)
454 >    call gather(q_group, q_group_Col, plan_group_col_3d)
455          
456      if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
457 <       call gather(u_l,u_l_Row,plan_row3d)
458 <       call gather(u_l,u_l_Col,plan_col3d)
457 >       call gather(u_l, u_l_Row, plan_atom_row_3d)
458 >       call gather(u_l, u_l_Col, plan_atom_col_3d)
459        
460 <       call gather(A,A_Row,plan_row_rotation)
461 <       call gather(A,A_Col,plan_col_rotation)
460 >       call gather(A, A_Row, plan_atom_row_rotation)
461 >       call gather(A, A_Col, plan_atom_col_rotation)
462      endif
463      
464   #endif
465 <
466 < !! Begin force loop timing:
465 >    
466 >    !! Begin force loop timing:
467   #ifdef PROFILE
468      call cpu_time(forceTimeInitial)
469      nloops = nloops + 1
470   #endif
471 <  
471 >    
472 >    loopEnd = PAIR_LOOP
473      if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
474 <       !! See if we need to update neighbor lists
475 <       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
476 <       !! if_mpi_gather_stuff_for_prepair
477 <       !! do_prepair_loop_if_needed
478 <       !! if_mpi_scatter_stuff_from_prepair
479 <       !! if_mpi_gather_stuff_from_prepair_to_main_loop
480 <    
481 < !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
474 >       loopStart = PREPAIR_LOOP
475 >    else
476 >       loopStart = PAIR_LOOP
477 >    endif
478 >
479 >    do loop = loopStart, loopEnd
480 >
481 >       ! See if we need to update neighbor lists
482 >       ! (but only on the first time through):
483 >       if (loop .eq. loopStart) then
484 >          call checkNeighborList(nGroups, q_group, listSkin, update_nlist)
485 >       endif
486 >      
487 >       if (update_nlist) then
488 >          !! save current configuration and construct neighbor list
489   #ifdef IS_MPI
490 <    
491 <    if (update_nlist) then
490 >          call saveNeighborList(nGroupsInRow, q_group)
491 > #else
492 >          call saveNeighborList(nGroups, q_group)
493 > #endif        
494 >          neighborListSize = size(list)
495 >          nlist = 0
496 >       endif
497        
498 <       !! save current configuration, construct neighbor list,
499 <       !! and calculate forces
500 <       call saveNeighborList(nlocal, q)
501 <      
502 <       neighborListSize = size(list)
503 <       nlist = 0      
504 <      
471 <       do i = 1, nrow
472 <          point(i) = nlist + 1
498 >       istart = 1
499 > #ifdef IS_MPI
500 >       iend = nGroupsInRow
501 > #else
502 >       iend = nGroups - 1
503 > #endif
504 >       outer: do i = istart, iend
505            
506 <          prepair_inner: do j = 1, ncol
507 <            
508 <             if (skipThisPair(i,j)) cycle prepair_inner
509 <            
510 <             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
511 <            
512 <             if (rijsq < rlistsq) then            
513 <                
514 <                nlist = nlist + 1
515 <                
516 <                if (nlist > neighborListSize) then
517 <                   call expandNeighborList(nlocal, listerror)
518 <                   if (listerror /= 0) then
519 <                      error = -1
520 <                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
521 <                      return
522 <                   end if
523 <                   neighborListSize = size(list)
524 <                endif
525 <                
526 <                list(nlist) = j
527 <                call do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot_local)                      
506 >          if (update_nlist) point(i) = nlist + 1
507 >          
508 >          n_in_i = groupStartRow(i+1) - groupStartRow(i)
509 >          
510 >          if (update_nlist) then
511 > #ifdef IS_MPI
512 >             jstart = 1
513 >             jend = nGroupsInCol
514 > #else
515 >             jstart = i+1
516 >             jend = nGroups
517 > #endif
518 >          else            
519 >             jstart = point(i)
520 >             jend = point(i+1) - 1
521 >             ! make sure group i has neighbors
522 >             if (jstart .gt. jend) cycle outer
523 >          endif
524 >          
525 >          do jnab = jstart, jend
526 >             if (update_nlist) then
527 >                j = jnab
528 >             else
529 >                j = list(jnab)
530               endif
497 –          enddo prepair_inner
498 –       enddo
531  
532 <       point(nrow + 1) = nlist + 1
533 <      
534 <    else  !! (of update_check)
503 <
504 <       ! use the list to find the neighbors
505 <       do i = 1, nrow
506 <          JBEG = POINT(i)
507 <          JEND = POINT(i+1) - 1
508 <          ! check thiat molecule i has neighbors
509 <          if (jbeg .le. jend) then
510 <            
511 <             do jnab = jbeg, jend
512 <                j = list(jnab)
513 <
514 <                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
515 <                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
516 <                     u_l, A, f, t, pot_local)
517 <
518 <             enddo
519 <          endif
520 <       enddo
521 <    endif
522 <    
532 > #ifdef IS_MPI
533 >             call get_interatomic_vector(q_group_Row(:,i), &
534 >                  q_group_Col(:,j), d_grp, rgrpsq)
535   #else
536 <    
537 <    if (update_nlist) then
538 <      
527 <       ! save current configuration, contruct neighbor list,
528 <       ! and calculate forces
529 <       call saveNeighborList(natoms, q)
530 <      
531 <       neighborListSize = size(list)
532 <  
533 <       nlist = 0
536 >             call get_interatomic_vector(q_group(:,i), &
537 >                  q_group(:,j), d_grp, rgrpsq)
538 > #endif
539  
540 <       do i = 1, natoms-1
541 <          point(i) = nlist + 1
542 <          
543 <          prepair_inner: do j = i+1, natoms
544 <            
545 <             if (skipThisPair(i,j))  cycle prepair_inner
546 <                          
547 <             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
548 <          
549 <
550 <             if (rijsq < rlistsq) then
551 <
552 <          
553 <                nlist = nlist + 1
554 <              
555 <                if (nlist > neighborListSize) then
556 <                   call expandNeighborList(natoms, listerror)
557 <                   if (listerror /= 0) then
558 <                      error = -1
554 <                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
555 <                      return
556 <                   end if
557 <                   neighborListSize = size(list)
540 >             if (rgrpsq < rlistsq) then
541 >                if (update_nlist) then
542 >                   nlist = nlist + 1
543 >                  
544 >                   if (nlist > neighborListSize) then
545 > #ifdef IS_MPI                
546 >                      call expandNeighborList(nGroupsInRow, listerror)
547 > #else
548 >                      call expandNeighborList(nGroups, listerror)
549 > #endif
550 >                      if (listerror /= 0) then
551 >                         error = -1
552 >                         write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
553 >                         return
554 >                      end if
555 >                      neighborListSize = size(list)
556 >                   endif
557 >                  
558 >                   list(nlist) = j
559                  endif
560                  
561 <                list(nlist) = j
561 >                if (loop .eq. PAIR_LOOP) then
562 >                   vij = 0.0d0
563 >                   fij(1:3) = 0.0d0
564 >                endif
565                  
566 <                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
567 <                        u_l, A, f, t, pot)
566 >                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
567 >                     in_switching_region)
568                  
569 <             endif
570 <          enddo prepair_inner
571 <       enddo
572 <      
573 <       point(natoms) = nlist + 1
574 <      
575 <    else !! (update)
576 <  
577 <       ! use the list to find the neighbors
578 <       do i = 1, natoms-1
579 <          JBEG = POINT(i)
580 <          JEND = POINT(i+1) - 1
581 <          ! check thiat molecule i has neighbors
582 <          if (jbeg .le. jend) then
583 <            
584 <             do jnab = jbeg, jend
581 <                j = list(jnab)
582 <
583 <                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
584 <                call do_prepair(i, j, rijsq, d, do_pot, do_stress, &
585 <                     u_l, A, f, t, pot)
586 <
587 <             enddo
588 <          endif
589 <       enddo
590 <    endif    
591 < #endif
592 <    !! Do rest of preforce calculations
593 <    !! do necessary preforce calculations  
594 <    call do_preforce(nlocal,pot)
595 <   ! we have already updated the neighbor list set it to false...
596 <   update_nlist = .false.
597 <    else
598 <       !! See if we need to update neighbor lists for non pre-pair
599 <       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
600 <    endif
601 <
602 <
603 <
604 <
605 <
606 < !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>>>>>>>>>>>>>>>>
607 <
608 <
609 <
610 <
611 <  
569 >                n_in_j = groupStartCol(j+1) - groupStartCol(j)
570 >                
571 >                do ia = groupStartRow(i), groupStartRow(i+1)-1
572 >                  
573 >                   atom1 = groupListRow(ia)
574 >                  
575 >                   inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
576 >                      
577 >                      atom2 = groupListCol(jb)
578 >                      
579 >                      if (skipThisPair(atom1, atom2)) cycle inner
580 >                      
581 >                      if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
582 >                         d_atm(1:3) = d_grp(1:3)
583 >                         ratmsq = rgrpsq
584 >                      else
585   #ifdef IS_MPI
586 <    
587 <    if (update_nlist) then
588 <       !! save current configuration, construct neighbor list,
589 <       !! and calculate forces
590 <       call saveNeighborList(nlocal, q)
591 <      
592 <       neighborListSize = size(list)
593 <       nlist = 0      
594 <      
595 <       do i = 1, nrow
596 <
597 <          point(i) = nlist + 1
598 <          
599 <          inner: do j = 1, ncol
600 <            
601 <             if (skipThisPair(i,j)) cycle inner
602 <            
603 <             call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
604 <            
605 <             if (rijsq < rlistsq) then            
586 >                         call get_interatomic_vector(q_Row(:,atom1), &
587 >                              q_Col(:,atom2), d_atm, ratmsq)
588 > #else
589 >                         call get_interatomic_vector(q(:,atom1), &
590 >                              q(:,atom2), d_atm, ratmsq)
591 > #endif
592 >                      endif
593 >                      if (loop .eq. PREPAIR_LOOP) then
594 > #ifdef IS_MPI                      
595 >                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
596 >                              rgrpsq, d_grp, do_pot, do_stress, &
597 >                              u_l, A, f, t, pot_local)
598 > #else
599 >                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
600 >                              rgrpsq, d_grp, do_pot, do_stress, &
601 >                              u_l, A, f, t, pot)
602 > #endif                                              
603 >                      else
604 > #ifdef IS_MPI                      
605 >                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
606 >                              do_pot, &
607 >                              u_l, A, f, t, pot_local, vpair, fpair)
608 > #else
609 >                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
610 >                              do_pot,  &
611 >                              u_l, A, f, t, pot, vpair, fpair)
612 > #endif
613 >                         vij = vij + vpair
614 >                         fij(1:3) = fij(1:3) + fpair(1:3)
615 >                      endif
616 >                   enddo inner
617 >                enddo
618                  
619 <                nlist = nlist + 1
620 <                
621 <                if (nlist > neighborListSize) then
622 <                   call expandNeighborList(nlocal, listerror)
623 <                   if (listerror /= 0) then
624 <                      error = -1
625 <                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
626 <                      return
627 <                   end if
628 <                   neighborListSize = size(list)
619 >                if (loop .eq. PAIR_LOOP) then
620 >                   if (in_switching_region) then
621 >                      swderiv = vij*dswdr/rgrp
622 >                      fij(1) = fij(1) + swderiv*d_grp(1)
623 >                      fij(2) = fij(2) + swderiv*d_grp(2)
624 >                      fij(3) = fij(3) + swderiv*d_grp(3)
625 >                      
626 >                      do ia=groupStartRow(i), groupStartRow(i+1)-1
627 >                         atom1=groupListRow(ia)
628 >                         mf = mfactRow(atom1)
629 > #ifdef IS_MPI
630 >                         f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
631 >                         f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
632 >                         f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
633 > #else
634 >                         f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
635 >                         f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
636 >                         f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
637 > #endif
638 >                      enddo
639 >                      
640 >                      do jb=groupStartCol(j), groupStartCol(j+1)-1
641 >                         atom2=groupListCol(jb)
642 >                         mf = mfactCol(atom2)
643 > #ifdef IS_MPI
644 >                         f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
645 >                         f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
646 >                         f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
647 > #else
648 >                         f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
649 >                         f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
650 >                         f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
651 > #endif
652 >                      enddo
653 >                   endif
654 >                  
655 >                   if (do_stress) call add_stress_tensor(d_grp, fij)
656                  endif
657 <                
658 <                list(nlist) = j
659 <                                
648 <                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
649 <                     u_l, A, f, t, pot_local)
650 <                
651 <             endif
652 <          enddo inner
653 <       enddo
654 <
655 <       point(nrow + 1) = nlist + 1
657 >             end if
658 >          enddo
659 >       enddo outer
660        
661 <    else  !! (of update_check)
662 <
663 <       ! use the list to find the neighbors
664 <       do i = 1, nrow
665 <          JBEG = POINT(i)
666 <          JEND = POINT(i+1) - 1
667 <          ! check thiat molecule i has neighbors
668 <          if (jbeg .le. jend) then
669 <            
670 <             do jnab = jbeg, jend
671 <                j = list(jnab)
668 <
669 <                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
670 <                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
671 <                     u_l, A, f, t, pot_local)
672 <
673 <             enddo
661 >       if (update_nlist) then
662 > #ifdef IS_MPI
663 >          point(nGroupsInRow + 1) = nlist + 1
664 > #else
665 >          point(nGroups) = nlist + 1
666 > #endif
667 >          if (loop .eq. PREPAIR_LOOP) then
668 >             ! we just did the neighbor list update on the first
669 >             ! pass, so we don't need to do it
670 >             ! again on the second pass
671 >             update_nlist = .false.                              
672            endif
673 <       enddo
676 <    endif
677 <    
678 < #else
679 <    
680 <    if (update_nlist) then
681 <
682 <       ! save current configuration, contruct neighbor list,
683 <       ! and calculate forces
684 <       call saveNeighborList(natoms, q)
685 <      
686 <       neighborListSize = size(list)
687 <  
688 <       nlist = 0
689 <      
690 <       do i = 1, natoms-1
691 <          point(i) = nlist + 1
692 <          
693 <          inner: do j = i+1, natoms
673 >       endif
674              
675 <             if (skipThisPair(i,j))  cycle inner
676 <                          
677 <             call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
698 <          
699 <
700 <             if (rijsq < rlistsq) then
701 <                
702 <                nlist = nlist + 1
703 <              
704 <                if (nlist > neighborListSize) then
705 <                   call expandNeighborList(natoms, listerror)
706 <                   if (listerror /= 0) then
707 <                      error = -1
708 <                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
709 <                      return
710 <                   end if
711 <                   neighborListSize = size(list)
712 <                endif
713 <                
714 <                list(nlist) = j
715 <                
716 <                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
717 <                        u_l, A, f, t, pot)
718 <                
719 <             endif
720 <          enddo inner
721 <       enddo
675 >       if (loop .eq. PREPAIR_LOOP) then
676 >          call do_preforce(nlocal, pot)
677 >       endif
678        
679 <       point(natoms) = nlist + 1
724 <      
725 <    else !! (update)
726 <      
727 <       ! use the list to find the neighbors
728 <       do i = 1, natoms-1
729 <          JBEG = POINT(i)
730 <          JEND = POINT(i+1) - 1
731 <          ! check thiat molecule i has neighbors
732 <          if (jbeg .le. jend) then
733 <            
734 <             do jnab = jbeg, jend
735 <                j = list(jnab)
736 <
737 <                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
738 <                call do_pair(i, j, rijsq, d, do_pot, do_stress, &
739 <                     u_l, A, f, t, pot)
740 <
741 <             enddo
742 <          endif
743 <       enddo
744 <    endif
679 >    enddo
680      
681 < #endif
747 <    
748 <    ! phew, done with main loop.
749 <
750 < !! Do timing
681 >    !! Do timing
682   #ifdef PROFILE
683      call cpu_time(forceTimeFinal)
684      forceTime = forceTime + forceTimeFinal - forceTimeInitial
685 < #endif
686 <
756 <
685 > #endif    
686 >    
687   #ifdef IS_MPI
688      !!distribute forces
689 <  
689 >    
690      f_temp = 0.0_dp
691 <    call scatter(f_Row,f_temp,plan_row3d)
691 >    call scatter(f_Row,f_temp,plan_atom_row_3d)
692      do i = 1,nlocal
693         f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
694      end do
695 <
695 >    
696      f_temp = 0.0_dp
697 <    call scatter(f_Col,f_temp,plan_col3d)
697 >    call scatter(f_Col,f_temp,plan_atom_col_3d)
698      do i = 1,nlocal
699         f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
700      end do
701      
702      if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
703         t_temp = 0.0_dp
704 <       call scatter(t_Row,t_temp,plan_row3d)
704 >       call scatter(t_Row,t_temp,plan_atom_row_3d)
705         do i = 1,nlocal
706            t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
707         end do
708         t_temp = 0.0_dp
709 <       call scatter(t_Col,t_temp,plan_col3d)
709 >       call scatter(t_Col,t_temp,plan_atom_col_3d)
710        
711         do i = 1,nlocal
712            t(1:3,i) = t(1:3,i) + t_temp(1:3,i)
# Line 785 | Line 715 | contains
715      
716      if (do_pot) then
717         ! scatter/gather pot_row into the members of my column
718 <       call scatter(pot_Row, pot_Temp, plan_row)
719 <
718 >       call scatter(pot_Row, pot_Temp, plan_atom_row)
719 >      
720         ! scatter/gather pot_local into all other procs
721         ! add resultant to get total pot
722         do i = 1, nlocal
# Line 794 | Line 724 | contains
724         enddo
725        
726         pot_Temp = 0.0_DP
727 <
728 <       call scatter(pot_Col, pot_Temp, plan_col)
727 >      
728 >       call scatter(pot_Col, pot_Temp, plan_atom_col)
729         do i = 1, nlocal
730            pot_local = pot_local + pot_Temp(i)
731         enddo
732 <
733 <    endif    
732 >      
733 >    endif
734   #endif
735 <
735 >    
736      if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
737        
738         if (FF_uses_RF .and. SIM_uses_RF) then
739 <
739 >          
740   #ifdef IS_MPI
741 <          call scatter(rf_Row,rf,plan_row3d)
742 <          call scatter(rf_Col,rf_Temp,plan_col3d)
741 >          call scatter(rf_Row,rf,plan_atom_row_3d)
742 >          call scatter(rf_Col,rf_Temp,plan_atom_col_3d)
743            do i = 1,nlocal
744               rf(1:3,i) = rf(1:3,i) + rf_Temp(1:3,i)
745            end do
746   #endif
747            
748            do i = 1, nLocal
749 <
749 >            
750               rfpot = 0.0_DP
751   #ifdef IS_MPI
752               me_i = atid_row(i)
753   #else
754               me_i = atid(i)
755   #endif
756 <
756 >            
757               if (PropertyMap(me_i)%is_DP) then
758 <
758 >                
759                  mu_i = PropertyMap(me_i)%dipole_moment
760 <
760 >                
761                  !! The reaction field needs to include a self contribution
762                  !! to the field:
763                  call accumulate_self_rf(i, mu_i, u_l)
# Line 844 | Line 774 | contains
774            enddo
775         endif
776      endif
777 <
778 <
777 >    
778 >    
779   #ifdef IS_MPI
780 <
780 >    
781      if (do_pot) then
782         pot = pot + pot_local
783         !! we assume the c code will do the allreduce to get the total potential
784         !! we could do it right here if we needed to...
785      endif
786 <
786 >    
787      if (do_stress) then
788 <      call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
788 >       call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
789              mpi_comm_world,mpi_err)
790         call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
791              mpi_comm_world,mpi_err)
792      endif
793 <
793 >    
794   #else
795 <
795 >    
796      if (do_stress) then
797         tau = tau_Temp
798         virial = virial_Temp
799      endif
800      
801   #endif
802 <    
873 <    
802 >      
803    end subroutine do_force_loop
804 +  
805 +  subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
806 +       u_l, A, f, t, pot, vpair, fpair)
807  
808 <  subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
809 <
810 <    real( kind = dp ) :: pot
808 >    real( kind = dp ) :: pot, vpair, sw
809 >    real( kind = dp ), dimension(3) :: fpair
810 >    real( kind = dp ), dimension(nLocal)   :: mfact
811      real( kind = dp ), dimension(3,nLocal) :: u_l
812 <    real (kind=dp), dimension(9,nLocal) :: A
813 <    real (kind=dp), dimension(3,nLocal) :: f
814 <    real (kind=dp), dimension(3,nLocal) :: t
812 >    real( kind = dp ), dimension(9,nLocal) :: A
813 >    real( kind = dp ), dimension(3,nLocal) :: f
814 >    real( kind = dp ), dimension(3,nLocal) :: t
815  
816 <    logical, intent(inout) :: do_pot, do_stress
816 >    logical, intent(inout) :: do_pot
817      integer, intent(in) :: i, j
818 <    real ( kind = dp ), intent(inout)    :: rijsq
818 >    real ( kind = dp ), intent(inout) :: rijsq
819      real ( kind = dp )                :: r
820      real ( kind = dp ), intent(inout) :: d(3)
821      integer :: me_i, me_j
822  
891 –
823      r = sqrt(rijsq)
824 +    vpair = 0.0d0
825 +    fpair(1:3) = 0.0d0
826  
827   #ifdef IS_MPI
828 <    if (tagRow(i) .eq. tagColumn(j)) then
829 <       write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
828 >    if (AtomRowToGlobal(i) .eq. AtomColToGlobal(j)) then
829 >       write(0,*) 'do_pair is doing', i , j, AtomRowToGlobal(i), AtomColToGlobal(j)
830      endif
831      me_i = atid_row(i)
832      me_j = atid_col(j)
# Line 905 | Line 838 | contains
838      if (FF_uses_LJ .and. SIM_uses_LJ) then
839        
840         if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) then
841 <          call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
841 >          !write(*,*) 'calling lj with'
842 >          !write(*,*) i, j, r, rijsq
843 >          !write(*,'(3es12.3)') d(1), d(2), d(3)
844 >          !write(*,'(3es12.3)') sw, vpair, pot
845 >          !write(*,*)
846 >
847 >          call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
848         endif
849        
850      endif
# Line 913 | Line 852 | contains
852      if (FF_uses_charges .and. SIM_uses_charges) then
853        
854         if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then
855 <          call do_charge_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
855 >          call do_charge_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
856         endif
857        
858      endif
# Line 921 | Line 860 | contains
860      if (FF_uses_dipoles .and. SIM_uses_dipoles) then
861        
862         if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
863 <          call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
864 <               do_pot, do_stress)
863 >          call do_dipole_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, u_l, f, t, &
864 >               do_pot)
865            if (FF_uses_RF .and. SIM_uses_RF) then
866 <             call accumulate_rf(i, j, r, u_l)
867 <             call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
866 >             call accumulate_rf(i, j, r, u_l, sw)
867 >             call rf_correct_forces(i, j, d, r, u_l, sw, f, fpair)
868            endif          
869         endif
870  
# Line 934 | Line 873 | contains
873      if (FF_uses_Sticky .and. SIM_uses_sticky) then
874  
875         if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
876 <          call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
877 <               do_pot, do_stress)
876 >          call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, A, f, t, &
877 >               do_pot)
878         endif
879  
880      endif
# Line 944 | Line 883 | contains
883      if (FF_uses_GB .and. SIM_uses_GB) then
884        
885         if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
886 <          call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
887 <               do_pot, do_stress)          
886 >          call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, u_l, f, t, &
887 >               do_pot)
888         endif
889  
890      endif
# Line 953 | Line 892 | contains
892      if (FF_uses_EAM .and. SIM_uses_EAM) then
893        
894         if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
895 <          call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
895 >          call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
896 >               do_pot)
897         endif
898        
899      endif
900 <
900 >    
901    end subroutine do_pair
902  
903 +  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
904 +       do_pot, do_stress, u_l, A, f, t, pot)
905  
906 <
965 <  subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
966 <   real( kind = dp ) :: pot
906 >   real( kind = dp ) :: pot, sw
907     real( kind = dp ), dimension(3,nLocal) :: u_l
908     real (kind=dp), dimension(9,nLocal) :: A
909     real (kind=dp), dimension(3,nLocal) :: f
# Line 971 | Line 911 | contains
911    
912     logical, intent(inout) :: do_pot, do_stress
913     integer, intent(in) :: i, j
914 <   real ( kind = dp ), intent(inout)    :: rijsq
915 <   real ( kind = dp )                :: r
916 <   real ( kind = dp ), intent(inout) :: d(3)
914 >   real ( kind = dp ), intent(inout)    :: rijsq, rcijsq
915 >   real ( kind = dp )                :: r, rc
916 >   real ( kind = dp ), intent(inout) :: d(3), dc(3)
917    
918     logical :: is_EAM_i, is_EAM_j
919    
920     integer :: me_i, me_j
921    
922 <   r = sqrt(rijsq)
922 >
923 >    r = sqrt(rijsq)
924 >    if (SIM_uses_molecular_cutoffs) then
925 >       rc = sqrt(rcijsq)
926 >    else
927 >       rc = r
928 >    endif
929    
930  
931   #ifdef IS_MPI
932 <   if (tagRow(i) .eq. tagColumn(j)) then
933 <      write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
932 >   if (AtomRowToGlobal(i) .eq. AtomColToGlobal(j)) then
933 >      write(0,*) 'do_prepair is doing', i , j, AtomRowToGlobal(i), AtomColToGlobal(j)
934     endif
935    
936     me_i = atid_row(i)
# Line 996 | Line 942 | contains
942     me_j = atid(j)
943    
944   #endif
945 <    
945 >  
946     if (FF_uses_EAM .and. SIM_uses_EAM) then
947 <
947 >      
948        if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
949             call calc_EAM_prepair_rho(i, j, d, r, rijsq )
950 <
950 >      
951     endif
952    
953   end subroutine do_prepair
954 <
1009 <
1010 <
1011 <
1012 <  subroutine do_preforce(nlocal,pot)
1013 <    integer :: nlocal
1014 <    real( kind = dp ) :: pot
1015 <
1016 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
1017 <       call calc_EAM_preforce_Frho(nlocal,pot)
1018 <    endif
1019 <
1020 <
1021 <  end subroutine do_preforce
1022 <  
1023 <  
1024 <  subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
1025 <    
1026 <    real (kind = dp), dimension(3) :: q_i
1027 <    real (kind = dp), dimension(3) :: q_j
1028 <    real ( kind = dp ), intent(out) :: r_sq
1029 <    real( kind = dp ) :: d(3), scaled(3)
1030 <    integer i
1031 <
1032 <    d(1:3) = q_j(1:3) - q_i(1:3)
1033 <
1034 <    ! Wrap back into periodic box if necessary
1035 <    if ( SIM_uses_PBC ) then
1036 <      
1037 <       if( .not.boxIsOrthorhombic ) then
1038 <          ! calc the scaled coordinates.
1039 <          
1040 <          scaled = matmul(HmatInv, d)
1041 <          
1042 <          ! wrap the scaled coordinates
1043 <
1044 <          scaled = scaled  - anint(scaled)
1045 <          
1046 <
1047 <          ! calc the wrapped real coordinates from the wrapped scaled
1048 <          ! coordinates
1049 <
1050 <          d = matmul(Hmat,scaled)
1051 <
1052 <       else
1053 <          ! calc the scaled coordinates.
1054 <          
1055 <          do i = 1, 3
1056 <             scaled(i) = d(i) * HmatInv(i,i)
1057 <            
1058 <             ! wrap the scaled coordinates
1059 <            
1060 <             scaled(i) = scaled(i) - anint(scaled(i))
1061 <            
1062 <             ! calc the wrapped real coordinates from the wrapped scaled
1063 <             ! coordinates
1064 <
1065 <             d(i) = scaled(i)*Hmat(i,i)
1066 <          enddo
1067 <       endif
1068 <      
1069 <    endif
1070 <    
1071 <    r_sq = dot_product(d,d)
1072 <    
1073 <  end subroutine get_interatomic_vector
1074 <  
1075 <  subroutine zero_work_arrays()
1076 <    
1077 < #ifdef IS_MPI
1078 <
1079 <    q_Row = 0.0_dp
1080 <    q_Col = 0.0_dp  
1081 <    
1082 <    u_l_Row = 0.0_dp
1083 <    u_l_Col = 0.0_dp
1084 <    
1085 <    A_Row = 0.0_dp
1086 <    A_Col = 0.0_dp
1087 <    
1088 <    f_Row = 0.0_dp
1089 <    f_Col = 0.0_dp
1090 <    f_Temp = 0.0_dp
1091 <      
1092 <    t_Row = 0.0_dp
1093 <    t_Col = 0.0_dp
1094 <    t_Temp = 0.0_dp
954 >
955  
956 <    pot_Row = 0.0_dp
957 <    pot_Col = 0.0_dp
958 <    pot_Temp = 0.0_dp
959 <
960 <    rf_Row = 0.0_dp
961 <    rf_Col = 0.0_dp
962 <    rf_Temp = 0.0_dp
963 <
964 < #endif
965 <
956 > subroutine do_preforce(nlocal,pot)
957 >   integer :: nlocal
958 >   real( kind = dp ) :: pot
959 >  
960 >   if (FF_uses_EAM .and. SIM_uses_EAM) then
961 >      call calc_EAM_preforce_Frho(nlocal,pot)
962 >   endif
963 >  
964 >  
965 > end subroutine do_preforce
966  
967 <    if (FF_uses_EAM .and. SIM_uses_EAM) then
968 <       call clean_EAM()
969 <    endif
970 <
971 <
972 <
973 <
974 <
975 <    rf = 0.0_dp
976 <    tau_Temp = 0.0_dp
977 <    virial_Temp = 0.0_dp
978 <  end subroutine zero_work_arrays
979 <  
980 <  function skipThisPair(atom1, atom2) result(skip_it)
981 <    integer, intent(in) :: atom1
982 <    integer, intent(in), optional :: atom2
983 <    logical :: skip_it
984 <    integer :: unique_id_1, unique_id_2
985 <    integer :: me_i,me_j
986 <    integer :: i
987 <
988 <    skip_it = .false.
989 <    
990 <    !! there are a number of reasons to skip a pair or a particle
991 <    !! mostly we do this to exclude atoms who are involved in short
992 <    !! range interactions (bonds, bends, torsions), but we also need
993 <    !! to exclude some overcounted interactions that result from
994 <    !! the parallel decomposition
995 <    
996 < #ifdef IS_MPI
997 <    !! in MPI, we have to look up the unique IDs for each atom
998 <    unique_id_1 = tagRow(atom1)
999 < #else
1000 <    !! in the normal loop, the atom numbers are unique
1001 <    unique_id_1 = atom1
1002 < #endif
1003 <
1004 <    !! We were called with only one atom, so just check the global exclude
1005 <    !! list for this atom
1006 <    if (.not. present(atom2)) then
1007 <       do i = 1, nExcludes_global
1008 <          if (excludesGlobal(i) == unique_id_1) then
1009 <             skip_it = .true.
1010 <             return
1011 <          end if
1012 <       end do
1013 <       return
1014 <    end if
1015 <    
1016 < #ifdef IS_MPI
1017 <    unique_id_2 = tagColumn(atom2)
1018 < #else
1019 <    unique_id_2 = atom2
1020 < #endif
967 >
968 > subroutine get_interatomic_vector(q_i, q_j, d, r_sq)
969 >  
970 >   real (kind = dp), dimension(3) :: q_i
971 >   real (kind = dp), dimension(3) :: q_j
972 >   real ( kind = dp ), intent(out) :: r_sq
973 >   real( kind = dp ) :: d(3), scaled(3)
974 >   integer i
975 >  
976 >   d(1:3) = q_j(1:3) - q_i(1:3)
977 >  
978 >   ! Wrap back into periodic box if necessary
979 >   if ( SIM_uses_PBC ) then
980 >      
981 >      if( .not.boxIsOrthorhombic ) then
982 >         ! calc the scaled coordinates.
983 >        
984 >         scaled = matmul(HmatInv, d)
985 >        
986 >         ! wrap the scaled coordinates
987 >        
988 >         scaled = scaled  - anint(scaled)
989 >        
990 >        
991 >         ! calc the wrapped real coordinates from the wrapped scaled
992 >         ! coordinates
993 >        
994 >         d = matmul(Hmat,scaled)
995 >        
996 >      else
997 >         ! calc the scaled coordinates.
998 >        
999 >         do i = 1, 3
1000 >            scaled(i) = d(i) * HmatInv(i,i)
1001 >            
1002 >            ! wrap the scaled coordinates
1003 >            
1004 >            scaled(i) = scaled(i) - anint(scaled(i))
1005 >            
1006 >            ! calc the wrapped real coordinates from the wrapped scaled
1007 >            ! coordinates
1008 >            
1009 >            d(i) = scaled(i)*Hmat(i,i)
1010 >         enddo
1011 >      endif
1012 >      
1013 >   endif
1014 >  
1015 >   r_sq = dot_product(d,d)
1016 >  
1017 > end subroutine get_interatomic_vector
1018 >
1019 > subroutine zero_work_arrays()
1020 >  
1021 > #ifdef IS_MPI
1022 >  
1023 >   q_Row = 0.0_dp
1024 >   q_Col = 0.0_dp
1025  
1026 +   q_group_Row = 0.0_dp
1027 +   q_group_Col = 0.0_dp  
1028 +  
1029 +   u_l_Row = 0.0_dp
1030 +   u_l_Col = 0.0_dp
1031 +  
1032 +   A_Row = 0.0_dp
1033 +   A_Col = 0.0_dp
1034 +  
1035 +   f_Row = 0.0_dp
1036 +   f_Col = 0.0_dp
1037 +   f_Temp = 0.0_dp
1038 +  
1039 +   t_Row = 0.0_dp
1040 +   t_Col = 0.0_dp
1041 +   t_Temp = 0.0_dp
1042 +  
1043 +   pot_Row = 0.0_dp
1044 +   pot_Col = 0.0_dp
1045 +   pot_Temp = 0.0_dp
1046 +  
1047 +   rf_Row = 0.0_dp
1048 +   rf_Col = 0.0_dp
1049 +   rf_Temp = 0.0_dp
1050 +  
1051 + #endif
1052 +
1053 +   if (FF_uses_EAM .and. SIM_uses_EAM) then
1054 +      call clean_EAM()
1055 +   endif
1056 +  
1057 +   rf = 0.0_dp
1058 +   tau_Temp = 0.0_dp
1059 +   virial_Temp = 0.0_dp
1060 + end subroutine zero_work_arrays
1061 +
1062 + function skipThisPair(atom1, atom2) result(skip_it)
1063 +   integer, intent(in) :: atom1
1064 +   integer, intent(in), optional :: atom2
1065 +   logical :: skip_it
1066 +   integer :: unique_id_1, unique_id_2
1067 +   integer :: me_i,me_j
1068 +   integer :: i
1069 +  
1070 +   skip_it = .false.
1071 +  
1072 +   !! there are a number of reasons to skip a pair or a particle
1073 +   !! mostly we do this to exclude atoms who are involved in short
1074 +   !! range interactions (bonds, bends, torsions), but we also need
1075 +   !! to exclude some overcounted interactions that result from
1076 +   !! the parallel decomposition
1077 +  
1078   #ifdef IS_MPI
1079 <    !! this situation should only arise in MPI simulations
1080 <    if (unique_id_1 == unique_id_2) then
1081 <       skip_it = .true.
1082 <       return
1083 <    end if
1168 <    
1169 <    !! this prevents us from doing the pair on multiple processors
1170 <    if (unique_id_1 < unique_id_2) then
1171 <       if (mod(unique_id_1 + unique_id_2,2) == 0) then
1172 <          skip_it = .true.
1173 <          return
1174 <       endif
1175 <    else                
1176 <       if (mod(unique_id_1 + unique_id_2,2) == 1) then
1177 <          skip_it = .true.
1178 <          return
1179 <       endif
1180 <    endif
1079 >   !! in MPI, we have to look up the unique IDs for each atom
1080 >   unique_id_1 = AtomRowToGlobal(atom1)
1081 > #else
1082 >   !! in the normal loop, the atom numbers are unique
1083 >   unique_id_1 = atom1
1084   #endif
1085 +  
1086 +   !! We were called with only one atom, so just check the global exclude
1087 +   !! list for this atom
1088 +   if (.not. present(atom2)) then
1089 +      do i = 1, nExcludes_global
1090 +         if (excludesGlobal(i) == unique_id_1) then
1091 +            skip_it = .true.
1092 +            return
1093 +         end if
1094 +      end do
1095 +      return
1096 +   end if
1097 +  
1098 + #ifdef IS_MPI
1099 +   unique_id_2 = AtomColToGlobal(atom2)
1100 + #else
1101 +   unique_id_2 = atom2
1102 + #endif
1103 +  
1104 + #ifdef IS_MPI
1105 +   !! this situation should only arise in MPI simulations
1106 +   if (unique_id_1 == unique_id_2) then
1107 +      skip_it = .true.
1108 +      return
1109 +   end if
1110 +  
1111 +   !! this prevents us from doing the pair on multiple processors
1112 +   if (unique_id_1 < unique_id_2) then
1113 +      if (mod(unique_id_1 + unique_id_2,2) == 0) then
1114 +         skip_it = .true.
1115 +         return
1116 +      endif
1117 +   else                
1118 +      if (mod(unique_id_1 + unique_id_2,2) == 1) then
1119 +         skip_it = .true.
1120 +         return
1121 +      endif
1122 +   endif
1123 + #endif
1124 +  
1125 +   !! the rest of these situations can happen in all simulations:
1126 +   do i = 1, nExcludes_global      
1127 +      if ((excludesGlobal(i) == unique_id_1) .or. &
1128 +           (excludesGlobal(i) == unique_id_2)) then
1129 +         skip_it = .true.
1130 +         return
1131 +      endif
1132 +   enddo
1133 +  
1134 +   do i = 1, nSkipsForAtom(unique_id_1)
1135 +      if (skipsForAtom(unique_id_1, i) .eq. unique_id_2) then
1136 +         skip_it = .true.
1137 +         return
1138 +      endif
1139 +   end do
1140 +  
1141 +   return
1142 + end function skipThisPair
1143  
1144 <    !! the rest of these situations can happen in all simulations:
1145 <    do i = 1, nExcludes_global      
1146 <       if ((excludesGlobal(i) == unique_id_1) .or. &
1147 <            (excludesGlobal(i) == unique_id_2)) then
1148 <          skip_it = .true.
1149 <          return
1150 <       endif
1151 <    enddo
1152 <
1153 <    do i = 1, nExcludes_local
1154 <       if (excludesLocal(1,i) == unique_id_1) then
1155 <          if (excludesLocal(2,i) == unique_id_2) then
1156 <             skip_it = .true.
1157 <             return
1158 <          endif
1159 <       else
1199 <          if (excludesLocal(1,i) == unique_id_2) then
1200 <             if (excludesLocal(2,i) == unique_id_1) then
1201 <                skip_it = .true.
1202 <                return
1203 <             endif
1204 <          endif
1205 <       endif
1206 <    end do
1207 <    
1208 <    return
1209 <  end function skipThisPair
1210 <
1211 <  function FF_UsesDirectionalAtoms() result(doesit)
1212 <    logical :: doesit
1213 <    doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
1214 <         FF_uses_GB .or. FF_uses_RF
1215 <  end function FF_UsesDirectionalAtoms
1216 <  
1217 <  function FF_RequiresPrepairCalc() result(doesit)
1218 <    logical :: doesit
1219 <    doesit = FF_uses_EAM
1220 <  end function FF_RequiresPrepairCalc
1221 <  
1222 <  function FF_RequiresPostpairCalc() result(doesit)
1223 <    logical :: doesit
1224 <    doesit = FF_uses_RF
1225 <  end function FF_RequiresPostpairCalc
1226 <  
1144 > function FF_UsesDirectionalAtoms() result(doesit)
1145 >   logical :: doesit
1146 >   doesit = FF_uses_dipoles .or. FF_uses_sticky .or. &
1147 >        FF_uses_GB .or. FF_uses_RF
1148 > end function FF_UsesDirectionalAtoms
1149 >
1150 > function FF_RequiresPrepairCalc() result(doesit)
1151 >   logical :: doesit
1152 >   doesit = FF_uses_EAM
1153 > end function FF_RequiresPrepairCalc
1154 >
1155 > function FF_RequiresPostpairCalc() result(doesit)
1156 >   logical :: doesit
1157 >   doesit = FF_uses_RF
1158 > end function FF_RequiresPostpairCalc
1159 >
1160   #ifdef PROFILE
1161 <  function getforcetime() result(totalforcetime)
1162 <    real(kind=dp) :: totalforcetime
1163 <    totalforcetime = forcetime
1164 <  end function getforcetime
1161 > function getforcetime() result(totalforcetime)
1162 >   real(kind=dp) :: totalforcetime
1163 >   totalforcetime = forcetime
1164 > end function getforcetime
1165   #endif
1166 +
1167 + !! This cleans componets of force arrays belonging only to fortran
1168  
1169 < !! This cleans componets of force arrays belonging only to fortran
1170 <
1169 > subroutine add_stress_tensor(dpair, fpair)
1170 >  
1171 >   real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1172 >  
1173 >   ! because the d vector is the rj - ri vector, and
1174 >   ! because fx, fy, fz are the force on atom i, we need a
1175 >   ! negative sign here:  
1176 >  
1177 >   tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1178 >   tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1179 >   tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1180 >   tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1181 >   tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1182 >   tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1183 >   tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1184 >   tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1185 >   tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1186 >  
1187 >   !write(*,'(6es12.3)')  fpair(1:3), tau_Temp(1), tau_Temp(5), tau_temp(9)
1188 >   virial_Temp = virial_Temp + &
1189 >        (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1190 >  
1191 > end subroutine add_stress_tensor
1192 >
1193   end module do_Forces
1194 +

Diff Legend

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