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 2394 by chrisfen, Sun Oct 23 21:08:08 2005 UTC vs.
Revision 2432 by chuckv, Tue Nov 15 16:01:06 2005 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.62 2005-10-23 21:08:02 chrisfen Exp $, $Date: 2005-10-23 21:08:02 $, $Name: not supported by cvs2svn $, $Revision: 1.62 $
48 > !! @version $Id: doForces.F90,v 1.68 2005-11-15 16:01:06 chuckv Exp $, $Date: 2005-11-15 16:01:06 $, $Name: not supported by cvs2svn $, $Revision: 1.68 $
49  
50  
51   module doForces
# Line 62 | Line 62 | module doForces
62    use shapes
63    use vector_class
64    use eam
65 +  use suttonchen
66    use status
67   #ifdef IS_MPI
68    use mpiSimulation
# Line 92 | Line 93 | module doForces
93    logical, save :: FF_uses_Dipoles
94    logical, save :: FF_uses_GayBerne
95    logical, save :: FF_uses_EAM
96 +  logical, save :: FF_uses_SC
97 +  logical, save :: FF_uses_MEAM
98 +
99  
100    logical, save :: SIM_uses_DirectionalAtoms
101    logical, save :: SIM_uses_EAM
102 +  logical, save :: SIM_uses_SC
103 +  logical, save :: SIM_uses_MEAM
104    logical, save :: SIM_requires_postpair_calc
105    logical, save :: SIM_requires_prepair_calc
106    logical, save :: SIM_uses_PBC
# Line 158 | Line 164 | contains
164      logical :: i_is_GB
165      logical :: i_is_EAM
166      logical :: i_is_Shape
167 +    logical :: i_is_SC
168 +    logical :: i_is_MEAM
169      logical :: j_is_LJ
170      logical :: j_is_Elect
171      logical :: j_is_Sticky
# Line 165 | Line 173 | contains
173      logical :: j_is_GB
174      logical :: j_is_EAM
175      logical :: j_is_Shape
176 +    logical :: j_is_SC
177 +    logical :: j_is_MEAM
178      real(kind=dp) :: myRcut
179  
180 +
181      status = 0  
182  
183      if (.not. associated(atypes)) then
# Line 204 | Line 215 | contains
215         call getElementProperty(atypes, i, "is_GayBerne", i_is_GB)
216         call getElementProperty(atypes, i, "is_EAM", i_is_EAM)
217         call getElementProperty(atypes, i, "is_Shape", i_is_Shape)
218 +       call getElementProperty(atypes, i, "is_SC", i_is_SC)
219 +       call getElementProperty(atypes, i, "is_MEAM", i_is_MEAM)
220  
221         do j = i, nAtypes
222  
# Line 217 | Line 230 | contains
230            call getElementProperty(atypes, j, "is_GayBerne", j_is_GB)
231            call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
232            call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
233 +          call getElementProperty(atypes, j, "is_SC", j_is_SC)
234 +          call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM)
235  
236            if (i_is_LJ .and. j_is_LJ) then
237               iHash = ior(iHash, LJ_PAIR)            
# Line 238 | Line 253 | contains
253               iHash = ior(iHash, EAM_PAIR)
254            endif
255  
256 +          if (i_is_SC .and. j_is_SC) then
257 +             iHash = ior(iHash, SC_PAIR)
258 +          endif
259 +
260            if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
261            if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
262            if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
# Line 569 | Line 588 | contains
588    subroutine setSimVariables()
589      SIM_uses_DirectionalAtoms = SimUsesDirectionalAtoms()
590      SIM_uses_EAM = SimUsesEAM()
591 +    SIM_uses_SC  = SimUsesSC()
592      SIM_requires_postpair_calc = SimRequiresPostpairCalc()
593      SIM_requires_prepair_calc = SimRequiresPrepairCalc()
594      SIM_uses_PBC = SimUsesPBC()
# Line 644 | Line 664 | contains
664      integer, intent(out) :: thisStat  
665      integer :: my_status, nMatches
666      integer, pointer :: MatchList(:) => null()
647 –    real(kind=dp) :: rcut, rrf, rt, dielect
667  
668      !! assume things are copacetic, unless they aren't
669      thisStat = 0
# Line 744 | Line 763 | contains
763      integer :: nlist
764      real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
765      real( kind = DP ) :: sw, dswdr, swderiv, mf
766 +    real( kind = DP ) :: rVal
767      real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
768      real(kind=dp) :: rfpot, mu_i, virial
769      integer :: me_i, me_j, n_in_i, n_in_j
# Line 754 | Line 774 | contains
774      integer :: propPack_i, propPack_j
775      integer :: loopStart, loopEnd, loop
776      integer :: iHash
777 <    integer :: ig
777 >    integer :: i1
778    
779  
780      !! initialize local variables  
# Line 902 | Line 922 | contains
922  
923                     list(nlist) = j
924                  endif
925 <
926 <                if (loop .eq. PAIR_LOOP) then
907 <                   vij = 0.0d0
908 <                   fij(1:3) = 0.0d0
909 <                endif
910 <
911 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
912 <                     in_switching_region)
913 <
914 <                n_in_j = groupStartCol(j+1) - groupStartCol(j)
915 <
916 <                do ia = groupStartRow(i), groupStartRow(i+1)-1
925 >                
926 >                if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
927  
928 <                   atom1 = groupListRow(ia)
929 <
930 <                   inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
931 <
932 <                      atom2 = groupListCol(jb)
933 <
934 <                      if (skipThisPair(atom1, atom2)) cycle inner
935 <
936 <                      if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
937 <                         d_atm(1:3) = d_grp(1:3)
938 <                         ratmsq = rgrpsq
939 <                      else
940 < #ifdef IS_MPI
941 <                         call get_interatomic_vector(q_Row(:,atom1), &
942 <                              q_Col(:,atom2), d_atm, ratmsq)
943 < #else
944 <                         call get_interatomic_vector(q(:,atom1), &
945 <                              q(:,atom2), d_atm, ratmsq)
946 < #endif
947 <                      endif
948 <
949 <                      if (loop .eq. PREPAIR_LOOP) then
950 < #ifdef IS_MPI                      
951 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
952 <                              rgrpsq, d_grp, do_pot, do_stress, &
953 <                              eFrame, A, f, t, pot_local)
954 < #else
945 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
946 <                              rgrpsq, d_grp, do_pot, do_stress, &
947 <                              eFrame, A, f, t, pot)
948 < #endif                                              
949 <                      else
950 < #ifdef IS_MPI                      
951 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
952 <                              do_pot, eFrame, A, f, t, pot_local, vpair, &
953 <                              fpair, d_grp, rgrp)
928 >                   if (loop .eq. PAIR_LOOP) then
929 >                      vij = 0.0d0
930 >                      fij(1:3) = 0.0d0
931 >                   endif
932 >                  
933 >                   call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
934 >                        in_switching_region)
935 >                  
936 >                   n_in_j = groupStartCol(j+1) - groupStartCol(j)
937 >                  
938 >                   do ia = groupStartRow(i), groupStartRow(i+1)-1
939 >                      
940 >                      atom1 = groupListRow(ia)
941 >                      
942 >                      inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
943 >                        
944 >                         atom2 = groupListCol(jb)
945 >                        
946 >                         if (skipThisPair(atom1, atom2))  cycle inner
947 >                        
948 >                         if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
949 >                            d_atm(1:3) = d_grp(1:3)
950 >                            ratmsq = rgrpsq
951 >                         else
952 > #ifdef IS_MPI
953 >                            call get_interatomic_vector(q_Row(:,atom1), &
954 >                                 q_Col(:,atom2), d_atm, ratmsq)
955   #else
956 <                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
957 <                              do_pot, eFrame, A, f, t, pot, vpair, fpair, &
957 <                              d_grp, rgrp)
956 >                            call get_interatomic_vector(q(:,atom1), &
957 >                                 q(:,atom2), d_atm, ratmsq)
958   #endif
959 +                         endif
960 +                        
961 +                         if (loop .eq. PREPAIR_LOOP) then
962 + #ifdef IS_MPI                      
963 +                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
964 +                                 rgrpsq, d_grp, do_pot, do_stress, &
965 +                                 eFrame, A, f, t, pot_local)
966 + #else
967 +                            call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
968 +                                 rgrpsq, d_grp, do_pot, do_stress, &
969 +                                 eFrame, A, f, t, pot)
970 + #endif                                              
971 +                         else
972 + #ifdef IS_MPI                      
973 +                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
974 +                                 do_pot, eFrame, A, f, t, pot_local, vpair, &
975 +                                 fpair, d_grp, rgrp)
976 + #else
977 +                            call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
978 +                                 do_pot, eFrame, A, f, t, pot, vpair, fpair, &
979 +                                 d_grp, rgrp)
980 + #endif
981 +                            vij = vij + vpair
982 +                            fij(1:3) = fij(1:3) + fpair(1:3)
983 +                         endif
984 +                      enddo inner
985 +                   enddo
986  
987 <                         vij = vij + vpair
988 <                         fij(1:3) = fij(1:3) + fpair(1:3)
989 <                      endif
990 <                   enddo inner
991 <                enddo
992 <
993 <                if (loop .eq. PAIR_LOOP) then
994 <                   if (in_switching_region) then
995 <                      swderiv = vij*dswdr/rgrp
996 <                      fij(1) = fij(1) + swderiv*d_grp(1)
997 <                      fij(2) = fij(2) + swderiv*d_grp(2)
998 <                      fij(3) = fij(3) + swderiv*d_grp(3)
999 <
1000 <                      do ia=groupStartRow(i), groupStartRow(i+1)-1
974 <                         atom1=groupListRow(ia)
975 <                         mf = mfactRow(atom1)
976 < #ifdef IS_MPI
977 <                         f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
978 <                         f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
979 <                         f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
987 >                   if (loop .eq. PAIR_LOOP) then
988 >                      if (in_switching_region) then
989 >                         swderiv = vij*dswdr/rgrp
990 >                         fij(1) = fij(1) + swderiv*d_grp(1)
991 >                         fij(2) = fij(2) + swderiv*d_grp(2)
992 >                         fij(3) = fij(3) + swderiv*d_grp(3)
993 >                        
994 >                         do ia=groupStartRow(i), groupStartRow(i+1)-1
995 >                            atom1=groupListRow(ia)
996 >                            mf = mfactRow(atom1)
997 > #ifdef IS_MPI
998 >                            f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
999 >                            f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
1000 >                            f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
1001   #else
1002 <                         f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1003 <                         f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1004 <                         f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1002 >                            f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
1003 >                            f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
1004 >                            f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
1005   #endif
1006 <                      enddo
1007 <
1008 <                      do jb=groupStartCol(j), groupStartCol(j+1)-1
1009 <                         atom2=groupListCol(jb)
1010 <                         mf = mfactCol(atom2)
1006 >                         enddo
1007 >                        
1008 >                         do jb=groupStartCol(j), groupStartCol(j+1)-1
1009 >                            atom2=groupListCol(jb)
1010 >                            mf = mfactCol(atom2)
1011   #ifdef IS_MPI
1012 <                         f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1013 <                         f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1014 <                         f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1012 >                            f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1013 >                            f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1014 >                            f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1015   #else
1016 <                         f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1017 <                         f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1018 <                         f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1016 >                            f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1017 >                            f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1018 >                            f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1019   #endif
1020 <                      enddo
1021 <                   endif
1020 >                         enddo
1021 >                      endif
1022  
1023 <                   if (do_stress) call add_stress_tensor(d_grp, fij)
1023 >                      if (do_stress) call add_stress_tensor(d_grp, fij)
1024 >                   endif
1025                  endif
1026 <             end if
1026 >             endif
1027            enddo
1028 <
1028 >          
1029         enddo outer
1030  
1031         if (update_nlist) then
# Line 1090 | Line 1112 | contains
1112            
1113            ! we loop only over the local atoms, so we don't need row and column
1114            ! lookups for the types
1115 <
1115 >          
1116            me_i = atid(i)
1117            
1118            ! is the atom electrostatic?  See if it would have an
1119            ! electrostatic interaction with itself
1120            iHash = InteractionHash(me_i,me_i)
1121 <          
1121 >
1122            if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1123   #ifdef IS_MPI
1124 <             call rf_self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1124 >             call self_self(i, eFrame, pot_local(ELECTROSTATIC_POT), &
1125                    t, do_pot)
1126   #else
1127 <             call rf_self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1127 >             call self_self(i, eFrame, pot(ELECTROSTATIC_POT), &
1128                    t, do_pot)
1129 + #endif
1130 +          endif
1131 +  
1132 +          
1133 +          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
1134 +            
1135 +             ! loop over the excludes to accumulate RF stuff we've
1136 +             ! left out of the normal pair loop
1137 +            
1138 +             do i1 = 1, nSkipsForAtom(i)
1139 +                j = skipsForAtom(i, i1)
1140 +                
1141 +                ! prevent overcounting of the skips
1142 +                if (i.lt.j) then
1143 +                   call get_interatomic_vector(q(:,i), &
1144 +                        q(:,j), d_atm, ratmsq)
1145 +                   rVal = dsqrt(ratmsq)
1146 +                   call get_switch(ratmsq, sw, dswdr, rVal, group_switch, &
1147 +                        in_switching_region)
1148 + #ifdef IS_MPI
1149 +                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1150 +                        vpair, pot_local(ELECTROSTATIC_POT), f, t, do_pot)
1151 + #else
1152 +                   call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
1153 +                        vpair, pot(ELECTROSTATIC_POT), f, t, do_pot)
1154   #endif
1155 +                endif
1156 +             enddo
1157            endif
1158         enddo
1159      endif
# Line 1170 | Line 1219 | contains
1219   #endif
1220  
1221      iHash = InteractionHash(me_i, me_j)
1222 <
1222 >    
1223      if ( iand(iHash, LJ_PAIR).ne.0 ) then
1224         call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1225              pot(VDW_POT), f, do_pot)
1226      endif
1227 <
1227 >    
1228      if ( iand(iHash, ELECTROSTATIC_PAIR).ne.0 ) then
1229         call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
1230              pot(ELECTROSTATIC_POT), eFrame, f, t, do_pot)
1231      endif
1232 <
1232 >    
1233      if ( iand(iHash, STICKY_PAIR).ne.0 ) then
1234         call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1235              pot(HB_POT), A, f, t, do_pot)
1236      endif
1237 <
1237 >    
1238      if ( iand(iHash, STICKYPOWER_PAIR).ne.0 ) then
1239         call do_sticky_power_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1240              pot(HB_POT), A, f, t, do_pot)
1241      endif
1242 <
1242 >    
1243      if ( iand(iHash, GAYBERNE_PAIR).ne.0 ) then
1244         call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1245              pot(VDW_POT), A, f, t, do_pot)
# Line 1200 | Line 1249 | contains
1249         call do_gb_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1250              pot(VDW_POT), A, f, t, do_pot)
1251      endif
1252 <
1252 >    
1253      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1254         call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1255              pot(METALLIC_POT), f, do_pot)
1256      endif
1257 <
1257 >    
1258      if ( iand(iHash, SHAPE_PAIR).ne.0 ) then      
1259         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1260              pot(VDW_POT), A, f, t, do_pot)
1261      endif
1262 <
1262 >    
1263      if ( iand(iHash, SHAPE_LJ).ne.0 ) then      
1264         call do_shape_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1265              pot(VDW_POT), A, f, t, do_pot)
1266      endif
1267 +
1268 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1269 +       call do_SC_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1270 +            pot(METALLIC_POT), f, do_pot)
1271 +    endif
1272 +
1273      
1274 +    
1275    end subroutine do_pair
1276  
1277    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
# Line 1251 | Line 1307 | contains
1307      if ( iand(iHash, EAM_PAIR).ne.0 ) then      
1308              call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1309      endif
1310 +
1311 +    if ( iand(iHash, SC_PAIR).ne.0 ) then      
1312 +            call calc_SC_prepair_rho(i, j, d, r, rijsq )
1313 +    endif
1314      
1315    end subroutine do_prepair
1316  
# Line 1262 | Line 1322 | contains
1322      if (FF_uses_EAM .and. SIM_uses_EAM) then
1323         call calc_EAM_preforce_Frho(nlocal,pot(METALLIC_POT))
1324      endif
1325 +    if (FF_uses_SC .and. SIM_uses_SC) then
1326 +       call calc_SC_preforce_Frho(nlocal,pot(METALLIC_POT))
1327 +    endif
1328  
1329  
1330    end subroutine do_preforce
# Line 1346 | Line 1409 | contains
1409      pot_Col = 0.0_dp
1410      pot_Temp = 0.0_dp
1411  
1349 –    rf_Row = 0.0_dp
1350 –    rf_Col = 0.0_dp
1351 –    rf_Temp = 0.0_dp
1352 –
1412   #endif
1413  
1414      if (FF_uses_EAM .and. SIM_uses_EAM) then
1415         call clean_EAM()
1416      endif
1417  
1359 –    rf = 0.0_dp
1418      tau_Temp = 0.0_dp
1419      virial_Temp = 0.0_dp
1420    end subroutine zero_work_arrays
# Line 1450 | Line 1508 | contains
1508  
1509    function FF_RequiresPrepairCalc() result(doesit)
1510      logical :: doesit
1511 <    doesit = FF_uses_EAM
1511 >    doesit = FF_uses_EAM .or. FF_uses_SC &
1512 >         .or. FF_uses_MEAM
1513    end function FF_RequiresPrepairCalc
1514  
1515   #ifdef PROFILE

Diff Legend

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