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 2787 by gezelter, Mon Jun 5 18:24:45 2006 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.83 2006-06-05 18:24:45 gezelter Exp $, $Date: 2006-06-05 18:24:45 $, $Name: not supported by cvs2svn $, $Revision: 1.83 $
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 89 | 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 105 | 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 117 | 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 147 | Line 151 | module doForces
151    end type gtypeCutoffs
152    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
153  
154 +  real(kind=dp), dimension(3) :: boxDipole
155  
156   contains
157  
# Line 571 | 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 609 | 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 636 | Line 655 | contains
655       SIM_requires_prepair_calc = SimRequiresPrepairCalc()
656       SIM_uses_PBC = SimUsesPBC()
657       SIM_uses_SC = SimUsesSC()
658 <    
658 >     SIM_uses_AtomicVirial = SimUsesAtomicVirial()
659 >
660       haveSIMvariables = .true.
661      
662       return
# Line 795 | 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 810 | Line 830 | contains
830      integer :: loopStart, loopEnd, loop
831      integer :: iHash
832      integer :: i1
833 <  
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  
# Line 1021 | 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) then
1092 +                               call add_stress_tensor(d_atm, fpair, tau)
1093 +                            endif
1094                           endif
1095                        enddo inner
1096                     enddo
# Line 1035 | 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
1063 –
1064 –                      if (do_stress) call add_stress_tensor(d_grp, fij)
1167                     endif
1168                  endif
1169               endif
# Line 1194 | 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
1201 –    
1314      if (do_pot) then
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,mpi_sum, &
1320 <            mpi_comm_world,mpi_err)            
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
1323 >        
1324 >    if (do_box_dipole) then
1325 >
1326   #ifdef SINGLE_PRECISION
1327 <       call mpi_allreduce(tau_Temp, tau, 9,mpi_real,mpi_sum, &
1328 <            mpi_comm_world,mpi_err)
1329 <       call mpi_allreduce(virial_Temp, virial,1,mpi_real,mpi_sum, &
1330 <            mpi_comm_world,mpi_err)
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 <       call mpi_allreduce(tau_Temp, tau, 9,mpi_double_precision,mpi_sum, &
1343 <            mpi_comm_world,mpi_err)
1344 <       call mpi_allreduce(virial_Temp, virial,1,mpi_double_precision,mpi_sum, &
1345 <            mpi_comm_world,mpi_err)
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 <    endif
1225 <    
1226 < #else
1227 <    
1228 <    if (do_stress) then
1229 <       tau = tau_Temp
1230 <       virial = virial_Temp
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 1479 | Line 1637 | contains
1637         call clean_EAM()
1638      endif
1639  
1482 –    tau_Temp = 0.0_dp
1483 –    virial_Temp = 0.0_dp
1640    end subroutine zero_work_arrays
1641  
1642    function skipThisPair(atom1, atom2) result(skip_it)
# Line 1585 | 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  
1606 –    virial_Temp = virial_Temp + &
1607 –         (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1608 –
1763    end subroutine add_stress_tensor
1764  
1765   end module doForces

Diff Legend

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