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 3126 by gezelter, Fri Apr 6 21:53:43 2007 UTC vs.
Revision 3397 by chuckv, Tue May 27 16:39:06 2008 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.85 2007-04-06 21:53:41 gezelter Exp $, $Date: 2007-04-06 21:53:41 $, $Name: not supported by cvs2svn $, $Revision: 1.85 $
48 > !! @version $Id: doForces.F90,v 1.93 2008-05-27 16:39:05 chuckv Exp $, $Date: 2008-05-27 16:39:05 $, $Name: not supported by cvs2svn $, $Revision: 1.93 $
49  
50  
51   module doForces
# Line 62 | Line 62 | module doForces
62    use shapes
63    use vector_class
64    use eam
65 +  use MetalNonMetal
66    use suttonchen
67    use status
68   #ifdef IS_MPI
# Line 96 | Line 97 | module doForces
97    logical, save :: FF_uses_GayBerne
98    logical, save :: FF_uses_EAM
99    logical, save :: FF_uses_SC
100 <  logical, save :: FF_uses_MEAM
100 >  logical, save :: FF_uses_MNM
101  
102  
103    logical, save :: SIM_uses_DirectionalAtoms
104    logical, save :: SIM_uses_EAM
105    logical, save :: SIM_uses_SC
106 <  logical, save :: SIM_uses_MEAM
106 >  logical, save :: SIM_uses_MNM
107    logical, save :: SIM_requires_postpair_calc
108    logical, save :: SIM_requires_prepair_calc
109    logical, save :: SIM_uses_PBC
# Line 113 | Line 114 | module doForces
114  
115    real(kind=dp), save :: defaultRcut, defaultRsw, largestRcut
116    real(kind=dp), save :: skinThickness
117 <  logical, save :: defaultDoShift
117 >  logical, save :: defaultDoShiftPot
118 >  logical, save :: defaultDoShiftFrc
119  
120    public :: init_FF
121    public :: setCutoffs
# Line 169 | Line 171 | contains
171      logical :: i_is_EAM
172      logical :: i_is_Shape
173      logical :: i_is_SC
172 –    logical :: i_is_MEAM
174      logical :: j_is_LJ
175      logical :: j_is_Elect
176      logical :: j_is_Sticky
# Line 178 | Line 179 | contains
179      logical :: j_is_EAM
180      logical :: j_is_Shape
181      logical :: j_is_SC
181 –    logical :: j_is_MEAM
182      real(kind=dp) :: myRcut
183  
184      if (.not. associated(atypes)) then
# Line 216 | Line 216 | contains
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)
219  
220         do j = i, nAtypes
221  
# Line 231 | Line 230 | contains
230            call getElementProperty(atypes, j, "is_EAM", j_is_EAM)
231            call getElementProperty(atypes, j, "is_Shape", j_is_Shape)
232            call getElementProperty(atypes, j, "is_SC", j_is_SC)
234 –          call getElementProperty(atypes, j, "is_MEAM", j_is_MEAM)
233  
234            if (i_is_LJ .and. j_is_LJ) then
235               iHash = ior(iHash, LJ_PAIR)            
# Line 260 | Line 258 | contains
258            if (i_is_GB .and. j_is_GB) iHash = ior(iHash, GAYBERNE_PAIR)
259            if (i_is_GB .and. j_is_LJ) iHash = ior(iHash, GAYBERNE_LJ)
260            if (i_is_LJ .and. j_is_GB) iHash = ior(iHash, GAYBERNE_LJ)
261 +        
262 +          if ((i_is_EAM.or.i_is_SC).and.(.not.(j_is_EAM.or.j_is_SC))) iHash = ior(iHash, MNM_PAIR)
263 +          if ((j_is_EAM.or.j_is_SC).and.(.not.(i_is_EAM.or.i_is_SC))) iHash = ior(iHash, MNM_PAIR)
264  
265            if (i_is_Shape .and. j_is_Shape) iHash = ior(iHash, SHAPE_PAIR)
266            if (i_is_Shape .and. j_is_LJ) iHash = ior(iHash, SHAPE_LJ)
# Line 562 | Line 563 | contains
563      haveGtypeCutoffMap = .true.
564     end subroutine createGtypeCutoffMap
565  
566 <   subroutine setCutoffs(defRcut, defRsw)
566 >   subroutine setCutoffs(defRcut, defRsw, defSP, defSF)
567  
568       real(kind=dp),intent(in) :: defRcut, defRsw
569 +     logical, intent(in) :: defSP, defSF
570       character(len = statusMsgSize) :: errMsg
571       integer :: localError
572  
573       defaultRcut = defRcut
574       defaultRsw = defRsw
575      
576 <     defaultDoShift = .false.
576 >     defaultDoShiftPot = defSP
577 >     defaultDoShiftFrc = defSF
578 >
579       if (abs(defaultRcut-defaultRsw) .lt. 0.0001) then
580 <        
581 <        write(errMsg, *) &
582 <             'cutoffRadius and switchingRadius are set to the same', newline &
583 <             // tab, 'value.  OOPSE will use shifted ', newline &
584 <             // tab, 'potentials instead of switching functions.'
585 <        
586 <        call handleInfo("setCutoffs", errMsg)
587 <        
588 <        defaultDoShift = .true.
589 <        
580 >        if (defaultDoShiftFrc) then
581 >           write(errMsg, *) &
582 >                'cutoffRadius and switchingRadius are set to the', newline &
583 >                // tab, 'same value.  OOPSE will use shifted force', newline &
584 >                // tab, 'potentials instead of switching functions.'
585 >          
586 >           call handleInfo("setCutoffs", errMsg)
587 >        else
588 >           write(errMsg, *) &
589 >                'cutoffRadius and switchingRadius are set to the', newline &
590 >                // tab, 'same value.  OOPSE will use shifted', newline &
591 >                // tab, 'potentials instead of switching functions.'
592 >          
593 >           call handleInfo("setCutoffs", errMsg)
594 >          
595 >           defaultDoShiftPot = .true.
596 >        endif
597 >                
598       endif
599      
600       localError = 0
601 <     call setLJDefaultCutoff( defaultRcut, defaultDoShift )
601 >     call setLJDefaultCutoff( defaultRcut, defaultDoShiftPot, &
602 >          defaultDoShiftFrc )
603       call setElectrostaticCutoffRadius( defaultRcut, defaultRsw )
604       call setCutoffEAM( defaultRcut )
605       call setCutoffSC( defaultRcut )
606 +     call setMnMDefaultCutoff( defaultRcut, defaultDoShiftPot, &
607 +          defaultDoShiftFrc )
608       call set_switch(defaultRsw, defaultRcut)
609       call setHmatDangerousRcutValue(defaultRcut)
610          
# Line 779 | Line 794 | contains
794  
795    !! Does force loop over i,j pairs. Calls do_pair to calculates forces.
796    !------------------------------------------------------------->
797 <  subroutine do_force_loop(q, q_group, A, eFrame, f, t, tau, pot, &
797 >  subroutine do_force_loop(q, q_group, A, eFrame, f, t, tau, pot, particle_pot,&
798         do_pot_c, do_stress_c, error)
799      !! Position array provided by C, dimensioned by getNlocal
800      real ( kind = dp ), dimension(3, nLocal) :: q
# Line 797 | Line 812 | contains
812      !! Stress Tensor
813      real( kind = dp), dimension(9) :: tau  
814      real ( kind = dp ),dimension(LR_POT_TYPES) :: pot
815 +    real( kind = dp ), dimension(nLocal) :: particle_pot
816      logical ( kind = 2) :: do_pot_c, do_stress_c
817      logical :: do_pot
818      logical :: do_stress
# Line 1079 | Line 1095 | contains
1095                              call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1096                                   do_pot, eFrame, A, f, t, pot_local, vpair, &
1097                                   fpair, d_grp, rgrp, rCut)
1098 +                            ! particle_pot will be accumulated from row & column
1099 +                            ! arrays later
1100   #else
1101                              call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
1102 <                                 do_pot, eFrame, A, f, t, pot, vpair, fpair, &
1103 <                                 d_grp, rgrp, rCut)
1102 >                                 do_pot, eFrame, A, f, t, pot, vpair, &
1103 >                                 fpair, d_grp, rgrp, rCut)
1104 >                            particle_pot(atom1) = particle_pot(atom1) + vpair*sw
1105 >                            particle_pot(atom2) = particle_pot(atom2) + vpair*sw
1106   #endif
1107                              vij = vij + vpair
1108                              fij(1) = fij(1) + fpair(1)
1109                              fij(2) = fij(2) + fpair(2)
1110                              fij(3) = fij(3) + fpair(3)
1111 <                            if (do_stress.and.SIM_uses_AtomicVirial) then
1111 >                            if (do_stress) then
1112                                 call add_stress_tensor(d_atm, fpair, tau)
1113                              endif
1114                           endif
# Line 1098 | Line 1118 | contains
1118                     if (loop .eq. PAIR_LOOP) then
1119                        if (in_switching_region) then
1120                           swderiv = vij*dswdr/rgrp
1121 <                         fij(1) = fij(1) + swderiv*d_grp(1)
1122 <                         fij(2) = fij(2) + swderiv*d_grp(2)
1123 <                         fij(3) = fij(3) + swderiv*d_grp(3)
1121 >                         fg = swderiv*d_grp
1122 >
1123 >                         fij(1) = fij(1) + fg(1)
1124 >                         fij(2) = fij(2) + fg(2)
1125 >                         fij(3) = fij(3) + fg(3)
1126                          
1127 +                         if (do_stress .and. (n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
1128 +                            call add_stress_tensor(d_atm, fg, tau)
1129 +                         endif  
1130 +                        
1131                           do ia=groupStartRow(i), groupStartRow(i+1)-1
1132                              atom1=groupListRow(ia)
1133                              mf = mfactRow(atom1)
# Line 1117 | Line 1143 | contains
1143                              f(2,atom1) = f(2,atom1) + fg(2)
1144                              f(3,atom1) = f(3,atom1) + fg(3)
1145   #endif
1146 <                            if (do_stress.and.SIM_uses_AtomicVirial) then
1147 <                               ! find the distance between the atom and the center of
1148 <                               ! the cutoff group:
1146 >                            if (n_in_i .gt. 1) then
1147 >                               if (do_stress.and.SIM_uses_AtomicVirial) then
1148 >                                  ! find the distance between the atom and the center of
1149 >                                  ! the cutoff group:
1150   #ifdef IS_MPI
1151 <                               call get_interatomic_vector(q_Row(:,atom1), &
1152 <                                    q_group_Row(:,i), dag, rag)
1151 >                                  call get_interatomic_vector(q_Row(:,atom1), &
1152 >                                       q_group_Row(:,i), dag, rag)
1153   #else
1154 <                               call get_interatomic_vector(q(:,atom1), &
1155 <                                    q_group(:,i), dag, rag)
1154 >                                  call get_interatomic_vector(q(:,atom1), &
1155 >                                       q_group(:,i), dag, rag)
1156   #endif
1157 <                               call add_stress_tensor(dag,fg,tau)
1157 >                                  call add_stress_tensor(dag,fg,tau)
1158 >                               endif
1159                              endif
1160                           enddo
1161                          
# Line 1146 | Line 1174 | contains
1174                              f(2,atom2) = f(2,atom2) + fg(2)
1175                              f(3,atom2) = f(3,atom2) + fg(3)
1176   #endif
1177 <                            if (do_stress.and.SIM_uses_AtomicVirial) then
1178 <                               ! find the distance between the atom and the center of
1179 <                               ! the cutoff group:
1177 >                            if (n_in_j .gt. 1) then
1178 >                               if (do_stress.and.SIM_uses_AtomicVirial) then
1179 >                                  ! find the distance between the atom and the center of
1180 >                                  ! the cutoff group:
1181   #ifdef IS_MPI
1182 <                               call get_interatomic_vector(q_Col(:,atom2), &
1183 <                                    q_group_Col(:,j), dag, rag)
1182 >                                  call get_interatomic_vector(q_Col(:,atom2), &
1183 >                                       q_group_Col(:,j), dag, rag)
1184   #else
1185 <                               call get_interatomic_vector(q(:,atom2), &
1186 <                                    q_group(:,j), dag, rag)
1185 >                                  call get_interatomic_vector(q(:,atom2), &
1186 >                                       q_group(:,j), dag, rag)
1187   #endif
1188 <                               call add_stress_tensor(dag,fg,tau)                              
1189 <                            endif
1190 <                            
1188 >                                  call add_stress_tensor(dag,fg,tau)                              
1189 >                               endif
1190 >                            endif                            
1191                           enddo
1192                        endif
1193 <                      if (do_stress.and.(.not.SIM_uses_AtomicVirial)) then
1194 <                         call add_stress_tensor(d_grp, fij, tau)
1195 <                      endif
1193 >                      !if (do_stress.and.(.not.SIM_uses_AtomicVirial)) then
1194 >                      !   call add_stress_tensor(d_grp, fij, tau)
1195 >                      !endif
1196                     endif
1197                  endif
1198               endif
# Line 1186 | Line 1215 | contains
1215         endif
1216  
1217         if (loop .eq. PREPAIR_LOOP) then
1218 + #ifdef IS_MPI
1219 +          call do_preforce(nlocal, pot_local)
1220 + #else
1221            call do_preforce(nlocal, pot)
1222 + #endif
1223         endif
1224  
1225      enddo
# Line 1238 | Line 1271 | contains
1271                 + pot_Temp(1:LR_POT_TYPES,i)
1272         enddo
1273  
1274 +       do i = 1,LR_POT_TYPES
1275 +          particle_pot(1:nlocal) = particle_pot(1:nlocal) + pot_Temp(i,1:nlocal)
1276 +       enddo
1277 +
1278         pot_Temp = 0.0_DP
1279 +
1280         do i = 1,LR_POT_TYPES
1281            call scatter(pot_Col(i,:), pot_Temp(i,:), plan_atom_col)
1282         end do
1283 +
1284         do i = 1, nlocal
1285            pot_local(1:LR_POT_TYPES) = pot_local(1:LR_POT_TYPES)&
1286                 + pot_Temp(1:LR_POT_TYPES,i)
1287         enddo
1288  
1289 +       do i = 1,LR_POT_TYPES
1290 +          particle_pot(1:nlocal) = particle_pot(1:nlocal) + pot_Temp(i,1:nlocal)
1291 +       enddo
1292 +
1293 +
1294      endif
1295   #endif
1296  
# Line 1397 | Line 1441 | contains
1441  
1442      real( kind = dp ) :: vpair, sw
1443      real( kind = dp ), dimension(LR_POT_TYPES) :: pot
1444 +    real( kind = dp ), dimension(nLocal) :: particle_pot
1445      real( kind = dp ), dimension(3) :: fpair
1446      real( kind = dp ), dimension(nLocal)   :: mfact
1447      real( kind = dp ), dimension(9,nLocal) :: eFrame
# Line 1411 | Line 1456 | contains
1456      real ( kind = dp ), intent(inout) :: d(3)
1457      real ( kind = dp ), intent(inout) :: d_grp(3)
1458      real ( kind = dp ), intent(inout) :: rCut
1459 <    real ( kind = dp ) :: r
1459 >    real ( kind = dp ) :: r, pair_pot
1460      real ( kind = dp ) :: a_k, b_k, c_k, d_k, dx
1461      integer :: me_i, me_j
1462      integer :: k
# Line 1483 | Line 1528 | contains
1528              pot(METALLIC_POT), f, do_pot)
1529      endif
1530      
1531 +    if ( iand(iHash, MNM_PAIR).ne.0 ) then      
1532 +       call do_mnm_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1533 +            pot(VDW_POT), A, f, t, do_pot)
1534 +    endif
1535 +
1536    end subroutine do_pair
1537  
1538    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, rCut, &
# Line 1728 | Line 1778 | contains
1778  
1779    function FF_RequiresPrepairCalc() result(doesit)
1780      logical :: doesit
1781 <    doesit = FF_uses_EAM .or. FF_uses_SC &
1732 <         .or. FF_uses_MEAM
1781 >    doesit = FF_uses_EAM .or. FF_uses_SC
1782    end function FF_RequiresPrepairCalc
1783  
1784   #ifdef PROFILE

Diff Legend

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