ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-3.0/src/UseTheForce/doForces.F90
(Generate patch)

Comparing trunk/OOPSE-3.0/src/UseTheForce/doForces.F90 (file contents):
Revision 1706 by gezelter, Thu Nov 4 16:20:28 2004 UTC vs.
Revision 2085 by gezelter, Tue Mar 8 21:05:46 2005 UTC

# Line 1 | Line 1
1 + !!
2 + !! Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 + !!
4 + !! The University of Notre Dame grants you ("Licensee") a
5 + !! non-exclusive, royalty free, license to use, modify and
6 + !! redistribute this software in source and binary code form, provided
7 + !! that the following conditions are met:
8 + !!
9 + !! 1. Acknowledgement of the program authors must be made in any
10 + !!    publication of scientific results based in part on use of the
11 + !!    program.  An acceptable form of acknowledgement is citation of
12 + !!    the article in which the program was described (Matthew
13 + !!    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 + !!    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 + !!    Parallel Simulation Engine for Molecular Dynamics,"
16 + !!    J. Comput. Chem. 26, pp. 252-271 (2005))
17 + !!
18 + !! 2. Redistributions of source code must retain the above copyright
19 + !!    notice, this list of conditions and the following disclaimer.
20 + !!
21 + !! 3. Redistributions in binary form must reproduce the above copyright
22 + !!    notice, this list of conditions and the following disclaimer in the
23 + !!    documentation and/or other materials provided with the
24 + !!    distribution.
25 + !!
26 + !! This software is provided "AS IS," without a warranty of any
27 + !! kind. All express or implied conditions, representations and
28 + !! warranties, including any implied warranty of merchantability,
29 + !! fitness for a particular purpose or non-infringement, are hereby
30 + !! excluded.  The University of Notre Dame and its licensors shall not
31 + !! be liable for any damages suffered by licensee as a result of
32 + !! using, modifying or distributing the software or its
33 + !! derivatives. In no event will the University of Notre Dame or its
34 + !! licensors be liable for any lost revenue, profit or data, or for
35 + !! direct, indirect, special, consequential, incidental or punitive
36 + !! damages, however caused and regardless of the theory of liability,
37 + !! arising out of the use of or inability to use software, even if the
38 + !! University of Notre Dame has been advised of the possibility of
39 + !! such damages.
40 + !!
41 +
42   !! doForces.F90
43   !! module doForces
44   !! Calculates Long Range forces.
45  
46   !! @author Charles F. Vardeman II
47   !! @author Matthew Meineke
48 < !! @version $Id: doForces.F90,v 1.7 2004-11-04 16:20:28 gezelter Exp $, $Date: 2004-11-04 16:20:28 $, $Name: not supported by cvs2svn $, $Revision: 1.7 $
48 > !! @version $Id: doForces.F90,v 1.11 2005-03-08 21:05:46 gezelter Exp $, $Date: 2005-03-08 21:05:46 $, $Name: not supported by cvs2svn $, $Revision: 1.11 $
49  
50 +
51   module doForces
52    use force_globals
53    use simulation
# Line 14 | Line 56 | module doForces
56    use switcheroo
57    use neighborLists  
58    use lj
59 <  use sticky_pair
60 <  use dipole_dipole
19 <  use charge_charge
59 >  use sticky
60 >  use electrostatic_module
61    use reaction_field
62    use gb_pair
63    use shapes
# Line 44 | Line 85 | module doForces
85    
86    logical, save :: FF_uses_DirectionalAtoms
87    logical, save :: FF_uses_LennardJones
88 <  logical, save :: FF_uses_Electrostatic
89 <  logical, save :: FF_uses_charges
90 <  logical, save :: FF_uses_dipoles
88 >  logical, save :: FF_uses_Electrostatics
89 >  logical, save :: FF_uses_Charges
90 >  logical, save :: FF_uses_Dipoles
91 >  logical, save :: FF_uses_Quadrupoles
92    logical, save :: FF_uses_sticky
93    logical, save :: FF_uses_GayBerne
94    logical, save :: FF_uses_EAM
# Line 59 | Line 101 | module doForces
101    logical, save :: SIM_uses_Electrostatics
102    logical, save :: SIM_uses_Charges
103    logical, save :: SIM_uses_Dipoles
104 +  logical, save :: SIM_uses_Quadrupoles
105    logical, save :: SIM_uses_Sticky
106    logical, save :: SIM_uses_GayBerne
107    logical, save :: SIM_uses_EAM
# Line 89 | Line 132 | module doForces
132       logical :: is_Electrostatic = .false.
133       logical :: is_Charge        = .false.
134       logical :: is_Dipole        = .false.
135 +     logical :: is_Quadrupole    = .false.
136       logical :: is_Sticky        = .false.
137       logical :: is_GayBerne      = .false.
138       logical :: is_EAM           = .false.
# Line 147 | Line 191 | contains
191         call getElementProperty(atypes, i, "is_Dipole", thisProperty)
192         PropertyMap(i)%is_Dipole = thisProperty
193  
194 +       call getElementProperty(atypes, i, "is_Quadrupole", thisProperty)
195 +       PropertyMap(i)%is_Quadrupole = thisProperty
196 +
197         call getElementProperty(atypes, i, "is_Sticky", thisProperty)
198         PropertyMap(i)%is_Sticky = thisProperty
199  
# Line 264 | Line 311 | contains
311    
312      FF_uses_DirectionalAtoms = .false.
313      FF_uses_LennardJones = .false.
314 <    FF_uses_Electrostatic = .false.
314 >    FF_uses_Electrostatics = .false.
315      FF_uses_Charges = .false.    
316      FF_uses_Dipoles = .false.
317      FF_uses_Sticky = .false.
# Line 284 | Line 331 | contains
331      call getMatchingElementList(atypes, "is_Electrostatic", .true., &
332           nMatches, MatchList)
333      if (nMatches .gt. 0) then
334 <       FF_uses_Electrostatic = .true.
334 >       FF_uses_Electrostatics = .true.
335      endif
336  
337      call getMatchingElementList(atypes, "is_Charge", .true., &
338           nMatches, MatchList)
339      if (nMatches .gt. 0) then
340 <       FF_uses_charges = .true.  
341 <       FF_uses_electrostatic = .true.
340 >       FF_uses_Charges = .true.  
341 >       FF_uses_Electrostatics = .true.
342      endif
343      
344      call getMatchingElementList(atypes, "is_Dipole", .true., &
345           nMatches, MatchList)
346      if (nMatches .gt. 0) then
347 <       FF_uses_dipoles = .true.
348 <       FF_uses_electrostatic = .true.
347 >       FF_uses_Dipoles = .true.
348 >       FF_uses_Electrostatics = .true.
349 >       FF_uses_DirectionalAtoms = .true.
350 >    endif
351 >
352 >    call getMatchingElementList(atypes, "is_Quadrupole", .true., &
353 >         nMatches, MatchList)
354 >    if (nMatches .gt. 0) then
355 >       FF_uses_Quadrupoles = .true.
356 >       FF_uses_Electrostatics = .true.
357         FF_uses_DirectionalAtoms = .true.
358      endif
359      
# Line 349 | Line 404 | contains
404         endif
405      endif
406  
407 <    if (FF_uses_sticky) then
408 <       call check_sticky_FF(my_status)
409 <       if (my_status /= 0) then
410 <          thisStat = -1
411 <          haveSaneForceField = .false.
412 <          return
413 <       end if
414 <    endif
407 >    !sticky module does not contain check_sticky_FF anymore
408 >    !if (FF_uses_sticky) then
409 >    !   call check_sticky_FF(my_status)
410 >    !   if (my_status /= 0) then
411 >    !      thisStat = -1
412 >    !      haveSaneForceField = .false.
413 >    !      return
414 >    !   end if
415 >    !endif
416  
417      if (FF_uses_EAM) then
418           call init_EAM_FF(my_status)
# Line 396 | Line 452 | contains
452  
453    !! Does force loop over i,j pairs. Calls do_pair to calculates forces.
454    !------------------------------------------------------------->
455 <  subroutine do_force_loop(q, q_group, A, u_l, f, t, tau, pot, &
455 >  subroutine do_force_loop(q, q_group, A, eFrame, f, t, tau, pot, &
456         do_pot_c, do_stress_c, error)
457      !! Position array provided by C, dimensioned by getNlocal
458      real ( kind = dp ), dimension(3, nLocal) :: q
# Line 405 | Line 461 | contains
461      !! Rotation Matrix for each long range particle in simulation.
462      real( kind = dp), dimension(9, nLocal) :: A    
463      !! Unit vectors for dipoles (lab frame)
464 <    real( kind = dp ), dimension(3,nLocal) :: u_l
464 >    real( kind = dp ), dimension(9,nLocal) :: eFrame
465      !! Force array provided by C, dimensioned by getNlocal
466      real ( kind = dp ), dimension(3,nLocal) :: f
467      !! Torsion array provided by C, dimensioned by getNlocal
# Line 480 | Line 536 | contains
536      call gather(q_group, q_group_Col, plan_group_col_3d)
537          
538      if (FF_UsesDirectionalAtoms() .and. SIM_uses_DirectionalAtoms) then
539 <       call gather(u_l, u_l_Row, plan_atom_row_3d)
540 <       call gather(u_l, u_l_Col, plan_atom_col_3d)
539 >       call gather(eFrame, eFrame_Row, plan_atom_row_rotation)
540 >       call gather(eFrame, eFrame_Col, plan_atom_col_rotation)
541        
542         call gather(A, A_Row, plan_atom_row_rotation)
543         call gather(A, A_Col, plan_atom_col_rotation)
# Line 627 | Line 683 | contains
683   #ifdef IS_MPI                      
684                           call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
685                                rgrpsq, d_grp, do_pot, do_stress, &
686 <                              u_l, A, f, t, pot_local)
686 >                              eFrame, A, f, t, pot_local)
687   #else
688                           call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
689                                rgrpsq, d_grp, do_pot, do_stress, &
690 <                              u_l, A, f, t, pot)
690 >                              eFrame, A, f, t, pot)
691   #endif                                              
692                        else
693   #ifdef IS_MPI                      
694                           call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
695                                do_pot, &
696 <                              u_l, A, f, t, pot_local, vpair, fpair)
696 >                              eFrame, A, f, t, pot_local, vpair, fpair)
697   #else
698                           call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
699                                do_pot,  &
700 <                              u_l, A, f, t, pot, vpair, fpair)
700 >                              eFrame, A, f, t, pot, vpair, fpair)
701   #endif
702  
703                           vij = vij + vpair
# Line 794 | Line 850 | contains
850                  
851                  !! The reaction field needs to include a self contribution
852                  !! to the field:
853 <                call accumulate_self_rf(i, mu_i, u_l)
853 >                call accumulate_self_rf(i, mu_i, eFrame)
854                  !! Get the reaction field contribution to the
855                  !! potential and torques:
856 <                call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
856 >                call reaction_field_final(i, mu_i, eFrame, rfpot, t, do_pot)
857   #ifdef IS_MPI
858                  pot_local = pot_local + rfpot
859   #else
# Line 837 | Line 893 | contains
893    end subroutine do_force_loop
894    
895    subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
896 <       u_l, A, f, t, pot, vpair, fpair)
896 >       eFrame, A, f, t, pot, vpair, fpair)
897  
898      real( kind = dp ) :: pot, vpair, sw
899      real( kind = dp ), dimension(3) :: fpair
900      real( kind = dp ), dimension(nLocal)   :: mfact
901 <    real( kind = dp ), dimension(3,nLocal) :: u_l
901 >    real( kind = dp ), dimension(9,nLocal) :: eFrame
902      real( kind = dp ), dimension(9,nLocal) :: A
903      real( kind = dp ), dimension(3,nLocal) :: f
904      real( kind = dp ), dimension(3,nLocal) :: t
# Line 866 | Line 922 | contains
922      me_j = atid(j)
923   #endif
924  
925 <    write(*,*) i, j, me_i, me_j
925 > !    write(*,*) i, j, me_i, me_j
926      
927      if (FF_uses_LennardJones .and. SIM_uses_LennardJones) then
928        
# Line 877 | Line 933 | contains
933        
934      endif
935      
936 <    if (FF_uses_charges .and. SIM_uses_charges) then
936 >    if (FF_uses_Electrostatics .and. SIM_uses_Electrostatics) then
937        
938 <       if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then
939 <          call do_charge_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
940 <               pot, f, do_pot)
938 >       if (PropertyMap(me_i)%is_Electrostatic .and. &
939 >            PropertyMap(me_j)%is_Electrostatic) then
940 >          call doElectrostaticPair(i, j, d, r, rijsq, sw, vpair, fpair, &
941 >               pot, eFrame, f, t, do_pot)
942         endif
943        
944 <    endif
945 <    
946 <    if (FF_uses_dipoles .and. SIM_uses_dipoles) then
947 <      
948 <       if ( PropertyMap(me_i)%is_Dipole .and. PropertyMap(me_j)%is_Dipole) then
949 <          call do_dipole_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
950 <               pot, u_l, f, t, do_pot)
894 <          if (FF_uses_RF .and. SIM_uses_RF) then
895 <             call accumulate_rf(i, j, r, u_l, sw)
896 <             call rf_correct_forces(i, j, d, r, u_l, sw, f, fpair)
944 >       if (FF_uses_dipoles .and. SIM_uses_dipoles) then      
945 >          if ( PropertyMap(me_i)%is_Dipole .and. &
946 >               PropertyMap(me_j)%is_Dipole) then
947 >             if (FF_uses_RF .and. SIM_uses_RF) then
948 >                call accumulate_rf(i, j, r, eFrame, sw)
949 >                call rf_correct_forces(i, j, d, r, eFrame, sw, f, fpair)
950 >             endif
951            endif
952         endif
899 –
953      endif
954  
955 +
956      if (FF_uses_Sticky .and. SIM_uses_sticky) then
957  
958         if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
# Line 914 | Line 968 | contains
968         if ( PropertyMap(me_i)%is_GayBerne .and. &
969              PropertyMap(me_j)%is_GayBerne) then
970            call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, &
971 <               pot, u_l, f, t, do_pot)
971 >               pot, A, f, t, do_pot)
972         endif
973        
974      endif
# Line 929 | Line 983 | contains
983      endif
984  
985  
986 <    write(*,*) PropertyMap(me_i)%is_Shape,PropertyMap(me_j)%is_Shape
986 > !    write(*,*) PropertyMap(me_i)%is_Shape,PropertyMap(me_j)%is_Shape
987  
988      if (FF_uses_Shapes .and. SIM_uses_Shapes) then
989         if ( PropertyMap(me_i)%is_Shape .and. &
# Line 943 | Line 997 | contains
997    end subroutine do_pair
998  
999    subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
1000 <       do_pot, do_stress, u_l, A, f, t, pot)
1000 >       do_pot, do_stress, eFrame, A, f, t, pot)
1001  
1002     real( kind = dp ) :: pot, sw
1003 <   real( kind = dp ), dimension(3,nLocal) :: u_l
1003 >   real( kind = dp ), dimension(9,nLocal) :: eFrame
1004     real (kind=dp), dimension(9,nLocal) :: A
1005     real (kind=dp), dimension(3,nLocal) :: f
1006     real (kind=dp), dimension(3,nLocal) :: t
# Line 1061 | Line 1115 | contains
1115     q_group_Row = 0.0_dp
1116     q_group_Col = 0.0_dp  
1117    
1118 <   u_l_Row = 0.0_dp
1119 <   u_l_Col = 0.0_dp
1118 >   eFrame_Row = 0.0_dp
1119 >   eFrame_Col = 0.0_dp
1120    
1121     A_Row = 0.0_dp
1122     A_Col = 0.0_dp
# Line 1179 | Line 1233 | contains
1233   function FF_UsesDirectionalAtoms() result(doesit)
1234     logical :: doesit
1235     doesit = FF_uses_DirectionalAtoms .or. FF_uses_Dipoles .or. &
1236 <        FF_uses_Sticky .or. FF_uses_GayBerne .or. FF_uses_Shapes
1236 >        FF_uses_Quadrupoles .or. FF_uses_Sticky .or. &
1237 >        FF_uses_GayBerne .or. FF_uses_Shapes
1238   end function FF_UsesDirectionalAtoms
1239  
1240   function FF_RequiresPrepairCalc() result(doesit)
# Line 1225 | Line 1280 | end module doForces
1280   end subroutine add_stress_tensor
1281  
1282   end module doForces
1228 –
1229 – !! Interfaces for C programs to module....
1230 –
1231 – subroutine initFortranFF(use_RF_c, thisStat)
1232 –    use doForces, ONLY: init_FF
1233 –    logical, intent(in) :: use_RF_c
1234 –
1235 –    integer, intent(out) :: thisStat  
1236 –    call init_FF(use_RF_c, thisStat)
1237 –
1238 – end subroutine initFortranFF
1239 –
1240 –  subroutine doForceloop(q, q_group, A, u_l, f, t, tau, pot, &
1241 –       do_pot_c, do_stress_c, error)
1242 –      
1243 –       use definitions, ONLY: dp
1244 –       use simulation
1245 –       use doForces, ONLY: do_force_loop
1246 –    !! Position array provided by C, dimensioned by getNlocal
1247 –    real ( kind = dp ), dimension(3, nLocal) :: q
1248 –    !! molecular center-of-mass position array
1249 –    real ( kind = dp ), dimension(3, nGroups) :: q_group
1250 –    !! Rotation Matrix for each long range particle in simulation.
1251 –    real( kind = dp), dimension(9, nLocal) :: A    
1252 –    !! Unit vectors for dipoles (lab frame)
1253 –    real( kind = dp ), dimension(3,nLocal) :: u_l
1254 –    !! Force array provided by C, dimensioned by getNlocal
1255 –    real ( kind = dp ), dimension(3,nLocal) :: f
1256 –    !! Torsion array provided by C, dimensioned by getNlocal
1257 –    real( kind = dp ), dimension(3,nLocal) :: t    
1258 –
1259 –    !! Stress Tensor
1260 –    real( kind = dp), dimension(9) :: tau  
1261 –    real ( kind = dp ) :: pot
1262 –    logical ( kind = 2) :: do_pot_c, do_stress_c
1263 –    integer :: error
1264 –    
1265 –    call do_force_loop(q, q_group, A, u_l, f, t, tau, pot, &
1266 –       do_pot_c, do_stress_c, error)
1267 –      
1268 – end subroutine doForceloop

Diff Legend

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