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 2722 by gezelter, Thu Apr 20 18:24:24 2006 UTC vs.
Revision 2917 by chrisfen, Mon Jul 3 13:18:43 2006 UTC

# Line 45 | Line 45
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.79 2006-04-20 18:24:24 gezelter Exp $, $Date: 2006-04-20 18:24:24 $, $Name: not supported by cvs2svn $, $Revision: 1.79 $
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 $
49  
50  
51   module doForces
# Line 64 | Line 64 | module doForces
64    use eam
65    use suttonchen
66    use status
67 –  use interpolation
67   #ifdef IS_MPI
68    use mpiSimulation
69   #endif
# Line 79 | Line 78 | module doForces
78  
79    INTEGER, PARAMETER:: PREPAIR_LOOP = 1
80    INTEGER, PARAMETER:: PAIR_LOOP    = 2
82 –  INTEGER, PARAMETER:: np = 500
81  
82    logical, save :: haveNeighborList = .false.
83    logical, save :: haveSIMvariables = .false.
# Line 91 | Line 89 | module doForces
89    logical, save :: haveElectrostaticSummationMethod = .false.
90    logical, save :: haveCutoffPolicy = .false.
91    logical, save :: VisitCutoffsAfterComputing = .false.
92 <  logical, save :: haveSplineSqrt = .false.
92 >  logical, save :: do_box_dipole = .false.
93  
94    logical, save :: FF_uses_DirectionalAtoms
95    logical, save :: FF_uses_Dipoles
# Line 120 | Line 118 | module doForces
118    public :: setCutoffs
119    public :: cWasLame
120    public :: setElectrostaticMethod
121 +  public :: setBoxDipole
122 +  public :: getBoxDipole
123    public :: setCutoffPolicy
124    public :: setSkinThickness
125    public :: do_force_loop
# Line 150 | Line 150 | module doForces
150    end type gtypeCutoffs
151    type(gtypeCutoffs), dimension(:,:), allocatable :: gtypeCutoffMap
152  
153 <  ! variables for the spline of the sqrt
154 <  type(cubicSpline), save :: splineSqrt
155 <  logical, save :: useSpline = .true.
156 <  
153 >  real(kind=dp), dimension(3) :: boxDipole
154  
155   contains
156  
# Line 430 | Line 427 | contains
427      !! largest cutoff for any atypes present in this group.  We also
428      !! create gtypes at this point.
429      
430 <    tol = 1.0d-6
430 >    tol = 1.0e-6_dp
431      nGroupTypesRow = 0
432      nGroupTypesCol = 0
433      do i = istart, iend      
# Line 594 | Line 591 | contains
591       call setCutoffSC( defaultRcut )
592       call set_switch(defaultRsw, defaultRcut)
593       call setHmatDangerousRcutValue(defaultRcut)
597 –     call setupSplineSqrt(defaultRcut)
594          
595       haveDefaultCutoffs = .true.
596       haveGtypeCutoffMap = .false.
# Line 618 | Line 614 | contains
614      
615     end subroutine setCutoffPolicy
616    
617 +   subroutine setBoxDipole()
618 +
619 +     do_box_dipole = .true.
620 +    
621 +   end subroutine setBoxDipole
622 +
623 +   subroutine getBoxDipole( box_dipole )
624 +
625 +     real(kind=dp), intent(inout), dimension(3) :: box_dipole
626 +
627 +     box_dipole = boxDipole
628 +
629 +   end subroutine getBoxDipole
630 +
631     subroutine setElectrostaticMethod( thisESM )
632  
633       integer, intent(in) :: thisESM
# Line 644 | Line 654 | contains
654       SIM_requires_prepair_calc = SimRequiresPrepairCalc()
655       SIM_uses_PBC = SimUsesPBC()
656       SIM_uses_SC = SimUsesSC()
657 <    
657 >
658       haveSIMvariables = .true.
659      
660       return
# Line 667 | Line 677 | contains
677      if (VisitCutoffsAfterComputing) then
678         call set_switch(largestRcut, largestRcut)      
679         call setHmatDangerousRcutValue(largestRcut)
670 –       call setLJsplineRmax(largestRcut)
680         call setCutoffEAM(largestRcut)
681         call setCutoffSC(largestRcut)
682         VisitCutoffsAfterComputing = .false.
683      endif
684  
676 –    if (.not. haveSplineSqrt) then
677 –       call setupSplineSqrt(largestRcut)
678 –    endif
679 –
685      if (.not. haveSIMvariables) then
686         call setSimVariables()
687      endif
# Line 823 | Line 828 | contains
828      integer :: loopStart, loopEnd, loop
829      integer :: iHash
830      integer :: i1
826 –  
831  
832 +    !! the variables for the box dipole moment
833 + #ifdef IS_MPI
834 +    integer :: pChgCount_local
835 +    integer :: nChgCount_local
836 +    real(kind=dp) :: pChg_local
837 +    real(kind=dp) :: nChg_local
838 +    real(kind=dp), dimension(3) :: pChgPos_local
839 +    real(kind=dp), dimension(3) :: nChgPos_local
840 +    real(kind=dp), dimension(3) :: dipVec_local
841 + #endif
842 +    integer :: pChgCount
843 +    integer :: nChgCount
844 +    real(kind=dp) :: pChg
845 +    real(kind=dp) :: nChg
846 +    real(kind=dp) :: chg_value
847 +    real(kind=dp), dimension(3) :: pChgPos
848 +    real(kind=dp), dimension(3) :: nChgPos
849 +    real(kind=dp), dimension(3) :: dipVec
850 +    real(kind=dp), dimension(3) :: chgVec
851 +
852 +    !! initialize box dipole variables
853 +    if (do_box_dipole) then
854 + #ifdef IS_MPI
855 +       pChg_local = 0.0_dp
856 +       nChg_local = 0.0_dp
857 +       pChgCount_local = 0
858 +       nChgCount_local = 0
859 +       do i=1, 3
860 +          pChgPos_local = 0.0_dp
861 +          nChgPos_local = 0.0_dp
862 +          dipVec_local = 0.0_dp
863 +       enddo
864 + #endif
865 +       pChg = 0.0_dp
866 +       nChg = 0.0_dp
867 +       pChgCount = 0
868 +       nChgCount = 0
869 +       chg_value = 0.0_dp
870 +      
871 +       do i=1, 3
872 +          pChgPos(i) = 0.0_dp
873 +          nChgPos(i) = 0.0_dp
874 +          dipVec(i) = 0.0_dp
875 +          chgVec(i) = 0.0_dp
876 +          boxDipole(i) = 0.0_dp
877 +       enddo
878 +    endif
879 +
880      !! initialize local variables  
881  
882   #ifdef IS_MPI
# Line 975 | Line 1027 | contains
1027  
1028                     rCut = gtypeCutoffMap(groupToGtypeRow(i),groupToGtypeCol(j))%rCut
1029                     if (loop .eq. PAIR_LOOP) then
1030 <                      vij = 0.0d0
1030 >                      vij = 0.0_dp
1031                        fij(1) = 0.0_dp
1032                        fij(2) = 0.0_dp
1033                        fij(3) = 0.0_dp
# Line 1195 | Line 1247 | contains
1247                  ! prevent overcounting of the skips
1248                  if (i.lt.j) then
1249                     call get_interatomic_vector(q(:,i), q(:,j), d_atm, ratmsq)
1250 <                   rVal = dsqrt(ratmsq)
1250 >                   rVal = sqrt(ratmsq)
1251                     call get_switch(ratmsq, sw, dswdr, rVal,in_switching_region)
1252   #ifdef IS_MPI
1253                     call rf_self_excludes(i, j, sw, eFrame, d_atm, rVal, &
# Line 1207 | Line 1259 | contains
1259                  endif
1260               enddo
1261            endif
1262 <       enddo
1263 <    endif
1212 <    
1262 >
1263 >          if (do_box_dipole) then
1264   #ifdef IS_MPI
1265 <    
1266 <    if (do_pot) then
1267 <       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision,mpi_sum, &
1265 >             call accumulate_box_dipole(i, eFrame, q(:,i), pChg_local, &
1266 >                  nChg_local, pChgPos_local, nChgPos_local, dipVec_local, &
1267 >                  pChgCount_local, nChgCount_local)
1268 > #else
1269 >             call accumulate_box_dipole(i, eFrame, q(:,i), pChg, nChg, &
1270 >                  pChgPos, nChgPos, dipVec, pChgCount, nChgCount)
1271 > #endif
1272 >          endif
1273 >       enddo
1274 >    endif
1275 >
1276 > #ifdef IS_MPI
1277 >    if (do_pot) then
1278 > #ifdef SINGLE_PRECISION
1279 >       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_real,mpi_sum, &
1280              mpi_comm_world,mpi_err)            
1281 + #else
1282 +       call mpi_allreduce(pot_local, pot, LR_POT_TYPES,mpi_double_precision, &
1283 +            mpi_sum, mpi_comm_world,mpi_err)            
1284 + #endif
1285      endif
1286      
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      
1301 +    if (do_box_dipole) then
1302 +
1303 + #ifdef SINGLE_PRECISION
1304 +       call mpi_allreduce(pChg_local, pChg, 1, mpi_real, mpi_sum, &
1305 +            mpi_comm_world, mpi_err)
1306 +       call mpi_allreduce(nChg_local, nChg, 1, mpi_real, mpi_sum, &
1307 +            mpi_comm_world, mpi_err)
1308 +       call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer, mpi_sum,&
1309 +            mpi_comm_world, mpi_err)
1310 +       call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer, mpi_sum,&
1311 +            mpi_comm_world, mpi_err)
1312 +       call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_real, mpi_sum, &
1313 +            mpi_comm_world, mpi_err)
1314 +       call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_real, mpi_sum, &
1315 +            mpi_comm_world, mpi_err)
1316 +       call mpi_allreduce(dipVec_local, dipVec, 3, mpi_real, mpi_sum, &
1317 +            mpi_comm_world, mpi_err)
1318   #else
1319 +       call mpi_allreduce(pChg_local, pChg, 1, mpi_double_precision, mpi_sum, &
1320 +            mpi_comm_world, mpi_err)
1321 +       call mpi_allreduce(nChg_local, nChg, 1, mpi_double_precision, mpi_sum, &
1322 +            mpi_comm_world, mpi_err)
1323 +       call mpi_allreduce(pChgCount_local, pChgCount, 1, mpi_integer,&
1324 +            mpi_sum, mpi_comm_world, mpi_err)
1325 +       call mpi_allreduce(nChgCount_local, nChgCount, 1, mpi_integer,&
1326 +            mpi_sum, mpi_comm_world, mpi_err)
1327 +       call mpi_allreduce(pChgPos_local, pChgPos, 3, mpi_double_precision, &
1328 +            mpi_sum, mpi_comm_world, mpi_err)
1329 +       call mpi_allreduce(nChgPos_local, nChgPos, 3, mpi_double_precision, &
1330 +            mpi_sum, mpi_comm_world, mpi_err)
1331 +       call mpi_allreduce(dipVec_local, dipVec, 3, mpi_double_precision, &
1332 +            mpi_sum, mpi_comm_world, mpi_err)
1333 + #endif
1334 +
1335 +    endif
1336 +
1337 + #else
1338      
1339      if (do_stress) then
1340         tau = tau_Temp
# Line 1232 | Line 1342 | contains
1342      endif
1343      
1344   #endif
1345 <    
1345 >
1346 >    if (do_box_dipole) then
1347 >       ! first load the accumulated dipole moment (if dipoles were present)
1348 >       boxDipole(1) = dipVec(1)
1349 >       boxDipole(2) = dipVec(2)
1350 >       boxDipole(3) = dipVec(3)
1351 >
1352 >       ! now include the dipole moment due to charges
1353 >       ! use the lesser of the positive and negative charge totals
1354 >       if (nChg .le. pChg) then
1355 >          chg_value = nChg
1356 >       else
1357 >          chg_value = pChg
1358 >       endif
1359 >      
1360 >       ! find the average positions
1361 >       if (pChgCount .gt. 0 .and. nChgCount .gt. 0) then
1362 >          pChgPos = pChgPos / pChgCount
1363 >          nChgPos = nChgPos / nChgCount
1364 >       endif
1365 >
1366 >       ! dipole is from the negative to the positive (physics notation)
1367 >       chgVec(1) = pChgPos(1) - nChgPos(1)
1368 >       chgVec(2) = pChgPos(2) - nChgPos(2)
1369 >       chgVec(3) = pChgPos(3) - nChgPos(3)
1370 >
1371 >       boxDipole(1) = boxDipole(1) + chgVec(1) * chg_value
1372 >       boxDipole(2) = boxDipole(2) + chgVec(2) * chg_value
1373 >       boxDipole(3) = boxDipole(3) + chgVec(3) * chg_value
1374 >
1375 >    endif
1376 >
1377    end subroutine do_force_loop
1378  
1379    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
# Line 1261 | Line 1402 | contains
1402  
1403      integer :: iHash
1404  
1405 <    if (useSpline) then
1406 <       call lookupUniformSpline(splineSqrt, rijsq, r)
1407 <    else
1408 <       r = sqrt(rijsq)
1268 <    endif
1405 >    r = sqrt(rijsq)
1406 >    
1407 >    vpair = 0.0_dp
1408 >    fpair(1:3) = 0.0_dp
1409  
1270 –    vpair = 0.0d0
1271 –    fpair(1:3) = 0.0d0
1272 –
1410   #ifdef IS_MPI
1411      me_i = atid_row(i)
1412      me_j = atid_col(j)
# Line 1306 | Line 1443 | contains
1443      endif
1444      
1445      if ( iand(iHash, GAYBERNE_LJ).ne.0 ) then
1446 <       call do_gb_lj_pair(i, j, d, r, rijsq, rcut, sw, vpair, fpair, &
1446 >       call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
1447              pot(VDW_POT), A, f, t, do_pot)
1448      endif
1449      
# Line 1350 | Line 1487 | contains
1487  
1488      integer :: me_i, me_j, iHash
1489  
1490 <    if (useSpline) then
1491 <       call lookupUniformSpline(splineSqrt, rijsq, r)
1355 <    else
1356 <       r = sqrt(rijsq)
1357 <    endif
1358 <
1490 >    r = sqrt(rijsq)
1491 >    
1492   #ifdef IS_MPI  
1493      me_i = atid_row(i)
1494      me_j = atid_col(j)  
# Line 1415 | Line 1548 | contains
1548            
1549            ! wrap the scaled coordinates
1550  
1551 <          scaled(1) = scaled(1) - dnint(scaled(1))
1552 <          scaled(2) = scaled(2) - dnint(scaled(2))
1553 <          scaled(3) = scaled(3) - dnint(scaled(3))
1551 >          scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1552 >          scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1553 >          scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1554  
1555            ! calc the wrapped real coordinates from the wrapped scaled
1556            ! coordinates
# Line 1435 | Line 1568 | contains
1568            
1569            ! wrap the scaled coordinates
1570            
1571 <          scaled(1) = scaled(1) - dnint(scaled(1))
1572 <          scaled(2) = scaled(2) - dnint(scaled(2))
1573 <          scaled(3) = scaled(3) - dnint(scaled(3))
1571 >          scaled(1) = scaled(1) - anint(scaled(1), kind=dp)
1572 >          scaled(2) = scaled(2) - anint(scaled(2), kind=dp)
1573 >          scaled(3) = scaled(3) - anint(scaled(3), kind=dp)
1574  
1575            ! calc the wrapped real coordinates from the wrapped scaled
1576            ! coordinates
# Line 1617 | Line 1750 | contains
1750  
1751    end subroutine add_stress_tensor
1752  
1620 –  subroutine setupSplineSqrt(rmax)
1621 –    real(kind=dp), intent(in) :: rmax
1622 –    real(kind=dp), dimension(np) :: xvals, yvals
1623 –    real(kind=dp) :: r2_1, r2_n, dx, r2
1624 –    integer :: i
1625 –
1626 –    r2_1 = 0.5d0
1627 –    r2_n = rmax*rmax
1628 –
1629 –    dx = (r2_n-r2_1) / dble(np-1)
1630 –    
1631 –    do i = 1, np
1632 –       r2 = r2_1 + dble(i-1)*dx
1633 –       xvals(i) = r2
1634 –       yvals(i) = dsqrt(r2)
1635 –    enddo
1636 –
1637 –    call newSpline(splineSqrt, xvals, yvals, .true.)
1638 –    
1639 –    haveSplineSqrt = .true.
1640 –    return
1641 –  end subroutine setupSplineSqrt
1642 –
1643 –  subroutine deleteSplineSqrt()
1644 –    call deleteSpline(splineSqrt)
1645 –    haveSplineSqrt = .false.
1646 –    return
1647 –  end subroutine deleteSplineSqrt
1648 –
1753   end module doForces

Diff Legend

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