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 2405 by chrisfen, Tue Nov 1 19:24:57 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.65 2005-11-01 19:24:54 chrisfen Exp $, $Date: 2005-11-01 19:24:54 $, $Name: not supported by cvs2svn $, $Revision: 1.65 $
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 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)
925 >                
926 >                if (rgrpsq < gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCutsq) then
927  
928 <                do ia = groupStartRow(i), groupStartRow(i+1)-1
929 <
930 <                   atom1 = groupListRow(ia)
931 <
932 <                   inner: do jb = groupStartCol(j), groupStartCol(j+1)-1
933 <
934 <                      atom2 = groupListCol(jb)
935 <    
936 <                      if (skipThisPair(atom1, atom2))  cycle inner
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 <                      if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
941 <                         d_atm(1:3) = d_grp(1:3)
942 <                         ratmsq = rgrpsq
943 <                      else
944 < #ifdef IS_MPI
945 <                         call get_interatomic_vector(q_Row(:,atom1), &
946 <                              q_Col(:,atom2), d_atm, ratmsq)
947 < #else
948 <                         call get_interatomic_vector(q(:,atom1), &
949 <                              q(:,atom2), d_atm, ratmsq)
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 get_interatomic_vector(q(:,atom1), &
957 >                                 q(:,atom2), d_atm, ratmsq)
958   #endif
959 <                      endif
960 <
961 <                      if (loop .eq. PREPAIR_LOOP) then
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)
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)
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
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)
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)
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 <
982 <                         vij = vij + vpair
983 <                         fij(1:3) = fij(1:3) + fpair(1:3)
984 <                      endif
985 <                   enddo inner
964 <                enddo
981 >                            vij = vij + vpair
982 >                            fij(1:3) = fij(1:3) + fpair(1:3)
983 >                         endif
984 >                      enddo inner
985 >                   enddo
986  
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)
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
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)
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
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
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
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
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 1108 | Line 1130 | contains
1130            endif
1131    
1132            
1133 < !!$          if (electrostaticSummationMethod.eq.REACTION_FIELD) then
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
# Line 1132 | Line 1154 | contains
1154   #endif
1155                  endif
1156               enddo
1157 < !!$          endif
1157 >          endif
1158         enddo
1159      endif
1160      
# Line 1242 | Line 1264 | contains
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  
# Line 1278 | 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 1289 | 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 1379 | Line 1415 | contains
1415         call clean_EAM()
1416      endif
1417  
1382 –    rf = 0.0_dp
1418      tau_Temp = 0.0_dp
1419      virial_Temp = 0.0_dp
1420    end subroutine zero_work_arrays
# Line 1473 | 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