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 2917 by chrisfen, Mon Jul 3 13:18:43 2006 UTC vs.
Revision 3126 by gezelter, Fri Apr 6 21:53:43 2007 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.84 2006-07-03 13:18:43 chrisfen Exp $, $Date: 2006-07-03 13:18:43 $, $Name: not supported by cvs2svn $, $Revision: 1.84 $
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 $
49  
50  
51   module doForces
# Line 106 | 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 613 | Line 614 | contains
614       haveGtypeCutoffMap = .false.
615      
616     end subroutine setCutoffPolicy
617 <  
617 >    
618     subroutine setBoxDipole()
619  
620       do_box_dipole = .true.
# Line 654 | Line 655 | contains
655       SIM_requires_prepair_calc = SimRequiresPrepairCalc()
656       SIM_uses_PBC = SimUsesPBC()
657       SIM_uses_SC = SimUsesSC()
658 +     SIM_uses_AtomicVirial = SimUsesAtomicVirial()
659  
660       haveSIMvariables = .true.
661      
# Line 813 | 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 1086 | Line 1088 | contains
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.and.SIM_uses_AtomicVirial) then
1092 +                               call add_stress_tensor(d_atm, fpair, tau)
1093 +                            endif
1094                           endif
1095                        enddo inner
1096                     enddo
# Line 1100 | 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 (do_stress.and.SIM_uses_AtomicVirial) then
1121 >                               ! find the distance between the atom and the center of
1122 >                               ! the cutoff group:
1123 > #ifdef IS_MPI
1124 >                               call get_interatomic_vector(q_Row(:,atom1), &
1125 >                                    q_group_Row(:,i), dag, rag)
1126 > #else
1127 >                               call get_interatomic_vector(q(:,atom1), &
1128 >                                    q_group(:,i), dag, rag)
1129   #endif
1130 +                               call add_stress_tensor(dag,fg,tau)
1131 +                            endif
1132                           enddo
1133                          
1134                           do jb=groupStartCol(j), groupStartCol(j+1)-1
1135                              atom2=groupListCol(jb)
1136                              mf = mfactCol(atom2)
1137 +                            ! fg is the force on atom jb due to cutoff group's
1138 +                            ! presence in switching region
1139 +                            fg = -swderiv*d_grp*mf
1140   #ifdef IS_MPI
1141 <                            f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
1142 <                            f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
1143 <                            f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
1141 >                            f_Col(1,atom2) = f_Col(1,atom2) + fg(1)
1142 >                            f_Col(2,atom2) = f_Col(2,atom2) + fg(2)
1143 >                            f_Col(3,atom2) = f_Col(3,atom2) + fg(3)
1144   #else
1145 <                            f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
1146 <                            f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
1147 <                            f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
1145 >                            f(1,atom2) = f(1,atom2) + fg(1)
1146 >                            f(2,atom2) = f(2,atom2) + fg(2)
1147 >                            f(3,atom2) = f(3,atom2) + fg(3)
1148 > #endif
1149 >                            if (do_stress.and.SIM_uses_AtomicVirial) then
1150 >                               ! find the distance between the atom and the center of
1151 >                               ! the cutoff group:
1152 > #ifdef IS_MPI
1153 >                               call get_interatomic_vector(q_Col(:,atom2), &
1154 >                                    q_group_Col(:,j), dag, rag)
1155 > #else
1156 >                               call get_interatomic_vector(q(:,atom2), &
1157 >                                    q_group(:,j), dag, rag)
1158   #endif
1159 +                               call add_stress_tensor(dag,fg,tau)                              
1160 +                            endif
1161 +                            
1162                           enddo
1163                        endif
1164 <
1165 <                      if (do_stress) call add_stress_tensor(d_grp, fij)
1164 >                      if (do_stress.and.(.not.SIM_uses_AtomicVirial)) then
1165 >                         call add_stress_tensor(d_grp, fij, tau)
1166 >                      endif
1167                     endif
1168                  endif
1169               endif
# Line 1283 | Line 1320 | contains
1320              mpi_sum, mpi_comm_world,mpi_err)            
1321   #endif
1322      endif
1323 <    
1287 <    if (do_stress) then
1288 < #ifdef SINGLE_PRECISION
1289 <       call mpi_allreduce(tau_Temp, tau, 9,mpi_real,mpi_sum, &
1290 <            mpi_comm_world,mpi_err)
1291 <       call mpi_allreduce(virial_Temp, virial,1,mpi_real,mpi_sum, &
1292 <            mpi_comm_world,mpi_err)
1293 < #else
1294 <       call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1295 <            mpi_comm_world,mpi_err)
1296 <       call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1297 <            mpi_comm_world,mpi_err)
1298 < #endif
1299 <    endif
1300 <    
1323 >        
1324      if (do_box_dipole) then
1325  
1326   #ifdef SINGLE_PRECISION
# Line 1333 | Line 1356 | contains
1356   #endif
1357  
1358      endif
1336 –
1337 – #else
1359      
1339 –    if (do_stress) then
1340 –       tau = tau_Temp
1341 –       virial = virial_Temp
1342 –    endif
1343 –    
1360   #endif
1361  
1362      if (do_box_dipole) then
# Line 1621 | Line 1637 | contains
1637         call clean_EAM()
1638      endif
1639  
1624 –    tau_Temp = 0.0_dp
1625 –    virial_Temp = 0.0_dp
1640    end subroutine zero_work_arrays
1641  
1642    function skipThisPair(atom1, atom2) result(skip_it)
# Line 1727 | 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  
1748 –    virial_Temp = virial_Temp + &
1749 –         (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1750 –
1763    end subroutine add_stress_tensor
1764  
1765   end module doForces

Diff Legend

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