ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/doForces.F90
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/doForces.F90 (file contents):
Revision 2512 by gezelter, Thu Dec 15 21:43:16 2005 UTC vs.
Revision 3127 by gezelter, Mon Apr 9 18:24:00 2007 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.71 2005-12-15 21:43:16 gezelter Exp $, $Date: 2005-12-15 21:43:16 $, $Name: not supported by cvs2svn $, $Revision: 1.71 $
48 > !! @version $Id: doForces.F90,v 1.86 2007-04-09 18:24:00 gezelter Exp $, $Date: 2007-04-09 18:24:00 $, $Name: not supported by cvs2svn $, $Revision: 1.86 $
49  
50  
51   module doForces
# Line 72 | Line 72 | module doForces
72    PRIVATE
73  
74   #define __FORTRAN90
75 – #include "UseTheForce/fSwitchingFunction.h"
75   #include "UseTheForce/fCutoffPolicy.h"
76   #include "UseTheForce/DarkSide/fInteractionMap.h"
77   #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
78  
80 –
79    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
80    INTEGER, PARAMETER:: PAIR_LOOP    = 2
81  
# Line 91 | Line 89 | module doForces
89    logical, save :: haveElectrostaticSummationMethod = .false.
90    logical, save :: haveCutoffPolicy = .false.
91    logical, save :: VisitCutoffsAfterComputing = .false.
92 +  logical, save :: do_box_dipole = .false.
93  
94    logical, save :: FF_uses_DirectionalAtoms
95    logical, save :: FF_uses_Dipoles
# Line 107 | Line 106 | module doForces
106    logical, save :: SIM_requires_postpair_calc
107    logical, save :: SIM_requires_prepair_calc
108    logical, save :: SIM_uses_PBC
109 +  logical, save :: SIM_uses_AtomicVirial
110  
111    integer, save :: electrostaticSummationMethod
112    integer, save :: cutoffPolicy = TRADITIONAL_CUTOFF_POLICY
# Line 119 | Line 119 | module doForces
119    public :: setCutoffs
120    public :: cWasLame
121    public :: setElectrostaticMethod
122 +  public :: setBoxDipole
123 +  public :: getBoxDipole
124    public :: setCutoffPolicy
125    public :: setSkinThickness
126    public :: do_force_loop
# Line 149 | Line 151 | contains
151    end type gtypeCutoffs
152    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
153  
154 +  real(kind=dp), dimension(3) :: boxDipole
155 +
156   contains
157  
158    subroutine createInteractionHash()
# Line 281 | Line 285 | contains
285      logical :: i_is_GB
286      logical :: i_is_EAM
287      logical :: i_is_Shape
288 +    logical :: i_is_SC
289      logical :: GtypeFound
290  
291      integer :: myStatus, nAtypes,  i, j, istart, iend, jstart, jend
# Line 310 | Line 315 | contains
315            call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
316            call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
317            call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
318 <          
318 >          call getElementProperty(atypes, i, "is_SC", i_is_SC)
319  
320            if (haveDefaultCutoffs) then
321               atypeMaxCutoff(i) = defaultRcut
# Line 343 | Line 348 | contains
348                  thisRcut = getShapeCut(i)
349                  if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
350               endif
351 +             if (i_is_SC) then
352 +                thisRcut = getSCCut(i)
353 +                if (thisRCut .gt. atypeMaxCutoff(i)) atypeMaxCutoff(i) = thisRCut
354 +             endif
355            endif
356                      
357            if (atypeMaxCutoff(i).gt.biggestAtypeCutoff) then
# Line 394 | Line 403 | contains
403         allocate(groupToGtypeCol(jend))
404      end if
405  
406 <    if(.not.associated(groupToGtypeCol)) then
407 <       allocate(groupToGtypeCol(jend))
406 >    if(.not.associated(groupMaxCutoffCol)) then
407 >       allocate(groupMaxCutoffCol(jend))
408      else
409 <       deallocate(groupToGtypeCol)
410 <       allocate(groupToGtypeCol(jend))
409 >       deallocate(groupMaxCutoffCol)
410 >       allocate(groupMaxCutoffCol(jend))
411      end if
412      if(.not.associated(gtypeMaxCutoffCol)) then
413         allocate(gtypeMaxCutoffCol(jend))
# Line 419 | Line 428 | contains
428      !! largest cutoff for any atypes present in this group.  We also
429      !! create gtypes at this point.
430      
431 <    tol = 1.0d-6
431 >    tol = 1.0e-6_dp
432      nGroupTypesRow = 0
433 <
433 >    nGroupTypesCol = 0
434      do i = istart, iend      
435         n_in_i = groupStartRow(i+1) - groupStartRow(i)
436         groupMaxCutoffRow(i) = 0.0_dp
# Line 567 | Line 576 | contains
576          
577          write(errMsg, *) &
578               'cutoffRadius and switchingRadius are set to the same', newline &
579 <             // tab, 'value.  OOPSE will use shifted ', newline &
579 >             // tab, 'value.  OOPSE will use shifted force van der Waals', newline &
580               // tab, 'potentials instead of switching functions.'
581          
582          call handleInfo("setCutoffs", errMsg)
# Line 575 | Line 584 | contains
584          defaultDoShift = .true.
585          
586       endif
587 <
587 >    
588       localError = 0
589       call setLJDefaultCutoff( defaultRcut, defaultDoShift )
590       call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
591 <     call setCutoffEAM( defaultRcut, localError)
592 <     if (localError /= 0) then
593 <       write(errMsg, *) 'An error has occured in setting the EAM cutoff'
594 <       call handleError("setCutoffs", errMsg)
595 <     end if
587 <     call set_switch(GROUP_SWITCH, defaultRsw, defaultRcut)
588 <
591 >     call setCutoffEAM( defaultRcut )
592 >     call setCutoffSC( defaultRcut )
593 >     call set_switch(defaultRsw, defaultRcut)
594 >     call setHmatDangerousRcutValue(defaultRcut)
595 >        
596       haveDefaultCutoffs = .true.
597       haveGtypeCutoffMap = .false.
598 +
599     end subroutine setCutoffs
600  
601     subroutine cWasLame()
# Line 606 | Line 614 | contains
614       haveGtypeCutoffMap = .false.
615      
616     end subroutine setCutoffPolicy
617 <  
617 >    
618 >   subroutine setBoxDipole()
619 >
620 >     do_box_dipole = .true.
621 >    
622 >   end subroutine setBoxDipole
623 >
624 >   subroutine getBoxDipole( box_dipole )
625 >
626 >     real(kind=dp), intent(inout), dimension(3) :: box_dipole
627 >
628 >     box_dipole = boxDipole
629 >
630 >   end subroutine getBoxDipole
631 >
632     subroutine setElectrostaticMethod( thisESM )
633  
634       integer, intent(in) :: thisESM
# Line 632 | Line 654 | contains
654       SIM_requires_postpair_calc = SimRequiresPostpairCalc()
655       SIM_requires_prepair_calc = SimRequiresPrepairCalc()
656       SIM_uses_PBC = SimUsesPBC()
657 <    
657 >     SIM_uses_SC = SimUsesSC()
658 >     SIM_uses_AtomicVirial = SimUsesAtomicVirial()
659 >
660       haveSIMvariables = .true.
661      
662       return
# Line 640 | Line 664 | contains
664  
665    subroutine doReadyCheck(error)
666      integer, intent(out) :: error
643 –
667      integer :: myStatus
668  
669      error = 0
# Line 653 | Line 676 | contains
676         call createGtypeCutoffMap()      
677      endif
678  
656 –
679      if (VisitCutoffsAfterComputing) then
680 <       call set_switch(GROUP_SWITCH, largestRcut, largestRcut)      
680 >       call set_switch(largestRcut, largestRcut)      
681 >       call setHmatDangerousRcutValue(largestRcut)
682 >       call setCutoffEAM(largestRcut)
683 >       call setCutoffSC(largestRcut)
684 >       VisitCutoffsAfterComputing = .false.
685      endif
686  
661 –
687      if (.not. haveSIMvariables) then
688         call setSimVariables()
689      endif
690  
666 –  !  if (.not. haveRlist) then
667 –  !     write(default_error, *) 'rList has not been set in doForces!'
668 –  !     error = -1
669 –  !     return
670 –  !  endif
671 –
691      if (.not. haveNeighborList) then
692         write(default_error, *) 'neighbor list has not been initialized in doForces!'
693         error = -1
694         return
695      end if
696 <
696 >    
697      if (.not. haveSaneForceField) then
698         write(default_error, *) 'Force Field is not sane in doForces!'
699         error = -1
700         return
701      end if
702 <
702 >    
703   #ifdef IS_MPI
704      if (.not. isMPISimSet()) then
705         write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
# Line 711 | Line 730 | contains
730      FF_uses_Dipoles = .false.
731      FF_uses_GayBerne = .false.
732      FF_uses_EAM = .false.
733 +    FF_uses_SC = .false.
734  
735      call getMatchingElementList(atypes, "is_Directional", .true., &
736           nMatches, MatchList)
# Line 727 | Line 747 | contains
747      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
748      if (nMatches .gt. 0) FF_uses_EAM = .true.
749  
750 +    call getMatchingElementList(atypes, "is_SC", .true., nMatches, MatchList)
751 +    if (nMatches .gt. 0) FF_uses_SC = .true.
752  
753 +
754      haveSaneForceField = .true.
755  
756      if (FF_uses_EAM) then
# Line 792 | Line 815 | contains
815      integer :: istart, iend
816      integer :: ia, jb, atom1, atom2
817      integer :: nlist
818 <    real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
818 >    real( kind = DP ) :: ratmsq, rgrpsq, rgrp, rag, vpair, vij
819      real( kind = DP ) :: sw, dswdr, swderiv, mf
820      real( kind = DP ) :: rVal
821 <    real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
822 <    real(kind=dp) :: rfpot, mu_i, virial
821 >    real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij, fg, dag
822 >    real(kind=dp) :: rfpot, mu_i
823      real(kind=dp):: rCut
824      integer :: me_i, me_j, n_in_i, n_in_j
825      logical :: is_dp_i
# Line 807 | Line 830 | contains
830      integer :: loopStart, loopEnd, loop
831      integer :: iHash
832      integer :: i1
810 –  
833  
834 +    !! the variables for the box dipole moment
835 + #ifdef IS_MPI
836 +    integer :: pChgCount_local
837 +    integer :: nChgCount_local
838 +    real(kind=dp) :: pChg_local
839 +    real(kind=dp) :: nChg_local
840 +    real(kind=dp), dimension(3) :: pChgPos_local
841 +    real(kind=dp), dimension(3) :: nChgPos_local
842 +    real(kind=dp), dimension(3) :: dipVec_local
843 + #endif
844 +    integer :: pChgCount
845 +    integer :: nChgCount
846 +    real(kind=dp) :: pChg
847 +    real(kind=dp) :: nChg
848 +    real(kind=dp) :: chg_value
849 +    real(kind=dp), dimension(3) :: pChgPos
850 +    real(kind=dp), dimension(3) :: nChgPos
851 +    real(kind=dp), dimension(3) :: dipVec
852 +    real(kind=dp), dimension(3) :: chgVec
853 +
854 +    !! initialize box dipole variables
855 +    if (do_box_dipole) then
856 + #ifdef IS_MPI
857 +       pChg_local = 0.0_dp
858 +       nChg_local = 0.0_dp
859 +       pChgCount_local = 0
860 +       nChgCount_local = 0
861 +       do i=1, 3
862 +          pChgPos_local = 0.0_dp
863 +          nChgPos_local = 0.0_dp
864 +          dipVec_local = 0.0_dp
865 +       enddo
866 + #endif
867 +       pChg = 0.0_dp
868 +       nChg = 0.0_dp
869 +       pChgCount = 0
870 +       nChgCount = 0
871 +       chg_value = 0.0_dp
872 +      
873 +       do i=1, 3
874 +          pChgPos(i) = 0.0_dp
875 +          nChgPos(i) = 0.0_dp
876 +          dipVec(i) = 0.0_dp
877 +          chgVec(i) = 0.0_dp
878 +          boxDipole(i) = 0.0_dp
879 +       enddo
880 +    endif
881 +
882      !! initialize local variables  
883  
884   #ifdef IS_MPI
# Line 954 | Line 1024 | contains
1024  
1025                     list(nlist) = j
1026                  endif
1027 <
958 <
959 <                
1027 >                
1028                  if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
1029  
1030                     rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
1031                     if (loop .eq. PAIR_LOOP) then
1032 <                      vij = 0.0d0
1033 <                      fij(1:3) = 0.0d0
1032 >                      vij = 0.0_dp
1033 >                      fij(1) = 0.0_dp
1034 >                      fij(2) = 0.0_dp
1035 >                      fij(3) = 0.0_dp
1036                     endif
1037                    
1038 <                   call get_switch(rgrpsq, sw, dswdr, rgrp, &
969 <                        group_switch, in_switching_region)
1038 >                   call get_switch(rgrpsq, sw, dswdr,rgrp, in_switching_region)
1039                    
1040                     n_in_j = groupStartCol(j+1) - groupStartCol(j)
1041                    
# Line 981 | Line 1050 | contains
1050                           if (skipThisPair(atom1, atom2))  cycle inner
1051                          
1052                           if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1053 <                            d_atm(1:3) = d_grp(1:3)
1053 >                            d_atm(1) = d_grp(1)
1054 >                            d_atm(2) = d_grp(2)
1055 >                            d_atm(3) = d_grp(3)
1056                              ratmsq = rgrpsq
1057                           else
1058   #ifdef IS_MPI
# Line 1014 | Line 1085 | contains
1085                                   d_grp, rgrp, rCut)
1086   #endif
1087                              vij = vij + vpair
1088 <                            fij(1:3) = fij(1:3) + fpair(1:3)
1088 >                            fij(1) = fij(1) + fpair(1)
1089 >                            fij(2) = fij(2) + fpair(2)
1090 >                            fij(3) = fij(3) + fpair(3)
1091 >                            if (do_stress) then
1092 >                               call add_stress_tensor(d_atm, fpair, tau)
1093 >                            endif
1094                           endif
1095                        enddo inner
1096                     enddo
# Line 1029 | Line 1105 | contains
1105                           do ia=groupStartRow(i), groupStartRow(i+1)-1
1106                              atom1=groupListRow(ia)
1107                              mf = mfactRow(atom1)
1108 +                            ! fg is the force on atom ia due to cutoff group's
1109 +                            ! presence in switching region
1110 +                            fg = swderiv*d_grp*mf
1111   #ifdef IS_MPI
1112 <                            f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
1113 <                            f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1114 <                            f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1112 >                            f_Row(1,atom1) = f_Row(1,atom1) + fg(1)
1113 >                            f_Row(2,atom1) = f_Row(2,atom1) + fg(2)
1114 >                            f_Row(3,atom1) = f_Row(3,atom1) + fg(3)
1115   #else
1116 <                            f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1117 <                            f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1118 <                            f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1116 >                            f(1,atom1) = f(1,atom1) + fg(1)
1117 >                            f(2,atom1) = f(2,atom1) + fg(2)
1118 >                            f(3,atom1) = f(3,atom1) + fg(3)
1119   #endif
1120 +                            if (n_in_i .gt. 1) then
1121 +                               if (do_stress.and.SIM_uses_AtomicVirial) then
1122 +                                  ! find the distance between the atom and the center of
1123 +                                  ! the cutoff group:
1124 + #ifdef IS_MPI
1125 +                                  call get_interatomic_vector(q_Row(:,atom1), &
1126 +                                       q_group_Row(:,i), dag, rag)
1127 + #else
1128 +                                  call get_interatomic_vector(q(:,atom1), &
1129 +                                       q_group(:,i), dag, rag)
1130 + #endif
1131 +                                  call add_stress_tensor(dag,fg,tau)
1132 +                               endif
1133 +                            endif
1134                           enddo
1135                          
1136                           do jb=groupStartCol(j), groupStartCol(j+1)-1
1137                              atom2=groupListCol(jb)
1138                              mf = mfactCol(atom2)
1139 +                            ! fg is the force on atom jb due to cutoff group's
1140 +                            ! presence in switching region
1141 +                            fg = -swderiv*d_grp*mf
1142   #ifdef IS_MPI
1143 <                            f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1144 <                            f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1145 <                            f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1143 >                            f_Col(1,atom2) = f_Col(1,atom2) + fg(1)
1144 >                            f_Col(2,atom2) = f_Col(2,atom2) + fg(2)
1145 >                            f_Col(3,atom2) = f_Col(3,atom2) + fg(3)
1146   #else
1147 <                            f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1148 <                            f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1149 <                            f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1147 >                            f(1,atom2) = f(1,atom2) + fg(1)
1148 >                            f(2,atom2) = f(2,atom2) + fg(2)
1149 >                            f(3,atom2) = f(3,atom2) + fg(3)
1150 > #endif
1151 >                            if (n_in_j .gt. 1) then
1152 >                               if (do_stress.and.SIM_uses_AtomicVirial) then
1153 >                                  ! find the distance between the atom and the center of
1154 >                                  ! the cutoff group:
1155 > #ifdef IS_MPI
1156 >                                  call get_interatomic_vector(q_Col(:,atom2), &
1157 >                                       q_group_Col(:,j), dag, rag)
1158 > #else
1159 >                                  call get_interatomic_vector(q(:,atom2), &
1160 >                                       q_group(:,j), dag, rag)
1161   #endif
1162 +                                  call add_stress_tensor(dag,fg,tau)                              
1163 +                               endif
1164 +                            endif                            
1165                           enddo
1166                        endif
1057 –
1058 –                      if (do_stress) call add_stress_tensor(d_grp, fij)
1167                     endif
1168                  endif
1169               endif
# Line 1175 | Line 1283 | contains
1283                  
1284                  ! prevent overcounting of the skips
1285                  if (i.lt.j) then
1286 <                   call get_interatomic_vector(q(:,i), &
1287 <                        q(:,j), d_atm, ratmsq)
1288 <                   rVal = dsqrt(ratmsq)
1181 <                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1182 <                        in_switching_region)
1286 >                   call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq)
1287 >                   rVal = sqrt(ratmsq)
1288 >                   call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region)
1289   #ifdef IS_MPI
1290                     call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1291                          vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
# Line 1190 | Line 1296 | contains
1296                  endif
1297               enddo
1298            endif
1299 +
1300 +          if (do_box_dipole) then
1301 + #ifdef IS_MPI
1302 +             call accumulate_box_dipole(i, eFrame, q(:,i), pChg_local, &
1303 +                  nChg_local, pChgPos_local, nChgPos_local, dipVec_local, &
1304 +                  pChgCount_local, nChgCount_local)
1305 + #else
1306 +             call accumulate_box_dipole(i, eFrame, q(:,i), pChg, nChg, &
1307 +                  pChgPos, nChgPos, dipVec, pChgCount, nChgCount)
1308 + #endif
1309 +          endif
1310         enddo
1311      endif
1312 <    
1312 >
1313   #ifdef IS_MPI
1197 –    
1314      if (do_pot) then
1315 <       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1315 > #ifdef SINGLE_PRECISION
1316 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_real,mpi_sum, &
1317              mpi_comm_world,mpi_err)            
1318 + #else
1319 +       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision, &
1320 +            mpi_sum, mpi_comm_world,mpi_err)            
1321 + #endif
1322      endif
1323 <    
1324 <    if (do_stress) then
1325 <       call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1326 <            mpi_comm_world,mpi_err)
1327 <       call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1328 <            mpi_comm_world,mpi_err)
1329 <    endif
1330 <    
1323 >        
1324 >    if (do_box_dipole) then
1325 >
1326 > #ifdef SINGLE_PRECISION
1327 >       call mpi_allreduce(pChg_local, pChg, 1, mpi_real, mpi_sum, &
1328 >            mpi_comm_world, mpi_err)
1329 >       call mpi_allreduce(nChg_local, nChg, 1, mpi_real, mpi_sum, &
1330 >            mpi_comm_world, mpi_err)
1331 >       call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer, mpi_sum,&
1332 >            mpi_comm_world, mpi_err)
1333 >       call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer, mpi_sum,&
1334 >            mpi_comm_world, mpi_err)
1335 >       call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_real, mpi_sum, &
1336 >            mpi_comm_world, mpi_err)
1337 >       call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_real, mpi_sum, &
1338 >            mpi_comm_world, mpi_err)
1339 >       call mpi_allreduce(dipVec_local, dipVec, 3, mpi_real, mpi_sum, &
1340 >            mpi_comm_world, mpi_err)
1341   #else
1342 <    
1343 <    if (do_stress) then
1344 <       tau = tau_Temp
1345 <       virial = virial_Temp
1342 >       call mpi_allreduce(pChg_local, pChg, 1, mpi_double_precision, mpi_sum, &
1343 >            mpi_comm_world, mpi_err)
1344 >       call mpi_allreduce(nChg_local, nChg, 1, mpi_double_precision, mpi_sum, &
1345 >            mpi_comm_world, mpi_err)
1346 >       call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer,&
1347 >            mpi_sum, mpi_comm_world, mpi_err)
1348 >       call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer,&
1349 >            mpi_sum, mpi_comm_world, mpi_err)
1350 >       call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_double_precision, &
1351 >            mpi_sum, mpi_comm_world, mpi_err)
1352 >       call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_double_precision, &
1353 >            mpi_sum, mpi_comm_world, mpi_err)
1354 >       call mpi_allreduce(dipVec_local, dipVec, 3, mpi_double_precision, &
1355 >            mpi_sum, mpi_comm_world, mpi_err)
1356 > #endif
1357 >
1358      endif
1359      
1360   #endif
1361 <    
1361 >
1362 >    if (do_box_dipole) then
1363 >       ! first load the accumulated dipole moment (if dipoles were present)
1364 >       boxDipole(1) = dipVec(1)
1365 >       boxDipole(2) = dipVec(2)
1366 >       boxDipole(3) = dipVec(3)
1367 >
1368 >       ! now include the dipole moment due to charges
1369 >       ! use the lesser of the positive and negative charge totals
1370 >       if (nChg .le. pChg) then
1371 >          chg_value = nChg
1372 >       else
1373 >          chg_value = pChg
1374 >       endif
1375 >      
1376 >       ! find the average positions
1377 >       if (pChgCount .gt. 0 .and. nChgCount .gt. 0) then
1378 >          pChgPos = pChgPos / pChgCount
1379 >          nChgPos = nChgPos / nChgCount
1380 >       endif
1381 >
1382 >       ! dipole is from the negative to the positive (physics notation)
1383 >       chgVec(1) = pChgPos(1) - nChgPos(1)
1384 >       chgVec(2) = pChgPos(2) - nChgPos(2)
1385 >       chgVec(3) = pChgPos(3) - nChgPos(3)
1386 >
1387 >       boxDipole(1) = boxDipole(1) + chgVec(1) * chg_value
1388 >       boxDipole(2) = boxDipole(2) + chgVec(2) * chg_value
1389 >       boxDipole(3) = boxDipole(3) + chgVec(3) * chg_value
1390 >
1391 >    endif
1392 >
1393    end subroutine do_force_loop
1394  
1395    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
# Line 1238 | Line 1412 | contains
1412      real ( kind = dp ), intent(inout) :: d_grp(3)
1413      real ( kind = dp ), intent(inout) :: rCut
1414      real ( kind = dp ) :: r
1415 +    real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1416      integer :: me_i, me_j
1417 +    integer :: k
1418  
1419      integer :: iHash
1420  
1421      r = sqrt(rijsq)
1422 <    vpair = 0.0d0
1423 <    fpair(1:3) = 0.0d0
1422 >    
1423 >    vpair = 0.0_dp
1424 >    fpair(1:3) = 0.0_dp
1425  
1426   #ifdef IS_MPI
1427      me_i = atid_row(i)
# Line 1282 | Line 1459 | contains
1459      endif
1460      
1461      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1462 <       call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1462 >       call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1463              pot(VDW_POT), A, f, t, do_pot)
1464      endif
1465      
# Line 1305 | Line 1482 | contains
1482         call do_SC_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1483              pot(METALLIC_POT), f, do_pot)
1484      endif
1308 –
1309 –    
1485      
1486    end subroutine do_pair
1487  
# Line 1329 | Line 1504 | contains
1504      integer :: me_i, me_j, iHash
1505  
1506      r = sqrt(rijsq)
1507 <
1507 >    
1508   #ifdef IS_MPI  
1509      me_i = atid_row(i)
1510      me_j = atid_col(j)  
# Line 1361 | Line 1536 | contains
1536      if (FF_uses_SC .and. SIM_uses_SC) then
1537         call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1538      endif
1364 –
1365 –
1539    end subroutine do_preforce
1540  
1541  
# Line 1374 | Line 1547 | contains
1547      real( kind = dp ) :: d(3), scaled(3)
1548      integer i
1549  
1550 <    d(1:3) = q_j(1:3) - q_i(1:3)
1550 >    d(1) = q_j(1) - q_i(1)
1551 >    d(2) = q_j(2) - q_i(2)
1552 >    d(3) = q_j(3) - q_i(3)
1553  
1554      ! Wrap back into periodic box if necessary
1555      if ( SIM_uses_PBC ) then
1556  
1557         if( .not.boxIsOrthorhombic ) then
1558            ! calc the scaled coordinates.
1559 +          ! scaled = matmul(HmatInv, d)
1560  
1561 <          scaled = matmul(HmatInv, d)
1562 <
1561 >          scaled(1) = HmatInv(1,1)*d(1) + HmatInv(1,2)*d(2) + HmatInv(1,3)*d(3)
1562 >          scaled(2) = HmatInv(2,1)*d(1) + HmatInv(2,2)*d(2) + HmatInv(2,3)*d(3)
1563 >          scaled(3) = HmatInv(3,1)*d(1) + HmatInv(3,2)*d(2) + HmatInv(3,3)*d(3)
1564 >          
1565            ! wrap the scaled coordinates
1566  
1567 <          scaled = scaled  - anint(scaled)
1567 >          scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1568 >          scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1569 >          scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1570  
1391 –
1571            ! calc the wrapped real coordinates from the wrapped scaled
1572            ! coordinates
1573 +          ! d = matmul(Hmat,scaled)
1574 +          d(1)= Hmat(1,1)*scaled(1) + Hmat(1,2)*scaled(2) + Hmat(1,3)*scaled(3)
1575 +          d(2)= Hmat(2,1)*scaled(1) + Hmat(2,2)*scaled(2) + Hmat(2,3)*scaled(3)
1576 +          d(3)= Hmat(3,1)*scaled(1) + Hmat(3,2)*scaled(2) + Hmat(3,3)*scaled(3)
1577  
1395 –          d = matmul(Hmat,scaled)
1396 –
1578         else
1579            ! calc the scaled coordinates.
1580  
1581 <          do i = 1, 3
1582 <             scaled(i) = d(i) * HmatInv(i,i)
1581 >          scaled(1) = d(1) * HmatInv(1,1)
1582 >          scaled(2) = d(2) * HmatInv(2,2)
1583 >          scaled(3) = d(3) * HmatInv(3,3)
1584 >          
1585 >          ! wrap the scaled coordinates
1586 >          
1587 >          scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1588 >          scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1589 >          scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1590  
1591 <             ! wrap the scaled coordinates
1591 >          ! calc the wrapped real coordinates from the wrapped scaled
1592 >          ! coordinates
1593  
1594 <             scaled(i) = scaled(i) - anint(scaled(i))
1594 >          d(1) = scaled(1)*Hmat(1,1)
1595 >          d(2) = scaled(2)*Hmat(2,2)
1596 >          d(3) = scaled(3)*Hmat(3,3)
1597  
1407 –             ! calc the wrapped real coordinates from the wrapped scaled
1408 –             ! coordinates
1409 –
1410 –             d(i) = scaled(i)*Hmat(i,i)
1411 –          enddo
1598         endif
1599  
1600      endif
1601  
1602 <    r_sq = dot_product(d,d)
1602 >    r_sq = d(1)*d(1) + d(2)*d(2) + d(3)*d(3)
1603  
1604    end subroutine get_interatomic_vector
1605  
# Line 1451 | Line 1637 | contains
1637         call clean_EAM()
1638      endif
1639  
1454 –    tau_Temp = 0.0_dp
1455 –    virial_Temp = 0.0_dp
1640    end subroutine zero_work_arrays
1641  
1642    function skipThisPair(atom1, atom2) result(skip_it)
# Line 1557 | Line 1741 | contains
1741  
1742    !! This cleans componets of force arrays belonging only to fortran
1743  
1744 <  subroutine add_stress_tensor(dpair, fpair)
1744 >  subroutine add_stress_tensor(dpair, fpair, tau)
1745  
1746      real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1747 +    real( kind = dp ), dimension(9), intent(inout) :: tau
1748  
1749      ! because the d vector is the rj - ri vector, and
1750      ! because fx, fy, fz are the force on atom i, we need a
1751      ! negative sign here:  
1752  
1753 <    tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1754 <    tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1755 <    tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1756 <    tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1757 <    tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1758 <    tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1759 <    tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1760 <    tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1761 <    tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1753 >    tau(1) = tau(1) - dpair(1) * fpair(1)
1754 >    tau(2) = tau(2) - dpair(1) * fpair(2)
1755 >    tau(3) = tau(3) - dpair(1) * fpair(3)
1756 >    tau(4) = tau(4) - dpair(2) * fpair(1)
1757 >    tau(5) = tau(5) - dpair(2) * fpair(2)
1758 >    tau(6) = tau(6) - dpair(2) * fpair(3)
1759 >    tau(7) = tau(7) - dpair(3) * fpair(1)
1760 >    tau(8) = tau(8) - dpair(3) * fpair(2)
1761 >    tau(9) = tau(9) - dpair(3) * fpair(3)
1762  
1578 –    virial_Temp = virial_Temp + &
1579 –         (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1580 –
1763    end subroutine add_stress_tensor
1764  
1765   end module doForces

Diff Legend

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