ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/do_Forces.F90
(Generate patch)

Comparing trunk/OOPSE/libmdtools/do_Forces.F90 (file contents):
Revision 669 by chuckv, Thu Aug 7 00:47:33 2003 UTC vs.
Revision 941 by gezelter, Tue Jan 13 23:01:43 2004 UTC

# Line 4 | Line 4
4  
5   !! @author Charles F. Vardeman II
6   !! @author Matthew Meineke
7 < !! @version $Id: do_Forces.F90,v 1.27 2003-08-07 00:47:33 chuckv Exp $, $Date: 2003-08-07 00:47:33 $, $Name: not supported by cvs2svn $, $Revision: 1.27 $
7 > !! @version $Id: do_Forces.F90,v 1.46 2004-01-13 23:01:43 gezelter Exp $, $Date: 2004-01-13 23:01:43 $, $Name: not supported by cvs2svn $, $Revision: 1.46 $
8  
9   module do_Forces
10    use force_globals
# Line 15 | Line 15 | module do_Forces
15    use lj
16    use sticky_pair
17    use dipole_dipole
18 +  use charge_charge
19    use reaction_field
20    use gb_pair
21    use vector_class
# Line 30 | Line 31 | module do_Forces
31   #define __FORTRAN90
32   #include "fForceField.h"
33  
34 <  logical, save :: do_forces_initialized = .false., haveRlist = .false.
34 >  logical, save :: haveRlist = .false.
35 >  logical, save :: haveNeighborList = .false.
36    logical, save :: havePolicies = .false.
37 +  logical, save :: haveSIMvariables = .false.
38 +  logical, save :: havePropertyMap = .false.
39 +  logical, save :: haveSaneForceField = .false.
40    logical, save :: FF_uses_LJ
41    logical, save :: FF_uses_sticky
42 +  logical, save :: FF_uses_charges
43    logical, save :: FF_uses_dipoles
44    logical, save :: FF_uses_RF
45    logical, save :: FF_uses_GB
46    logical, save :: FF_uses_EAM
47 +  logical, save :: SIM_uses_LJ
48 +  logical, save :: SIM_uses_sticky
49 +  logical, save :: SIM_uses_charges
50 +  logical, save :: SIM_uses_dipoles
51 +  logical, save :: SIM_uses_RF
52 +  logical, save :: SIM_uses_GB
53 +  logical, save :: SIM_uses_EAM
54 +  logical, save :: SIM_requires_postpair_calc
55 +  logical, save :: SIM_requires_prepair_calc
56 +  logical, save :: SIM_uses_directional_atoms
57 +  logical, save :: SIM_uses_PBC
58  
59    real(kind=dp), save :: rlist, rlistsq
60  
# Line 45 | Line 62 | contains
62    public :: do_force_loop
63    public :: setRlistDF
64  
65 +
66 + #ifdef PROFILE
67 +  public :: getforcetime
68 +  real, save :: forceTime = 0
69 +  real :: forceTimeInitial, forceTimeFinal
70 +  integer :: nLoops
71 + #endif
72 +
73 +  type :: Properties
74 +     logical :: is_lj     = .false.
75 +     logical :: is_sticky = .false.
76 +     logical :: is_dp     = .false.
77 +     logical :: is_gb     = .false.
78 +     logical :: is_eam    = .false.
79 +     logical :: is_charge = .false.
80 +     real(kind=DP) :: charge = 0.0_DP
81 +     real(kind=DP) :: dipole_moment = 0.0_DP
82 +  end type Properties
83 +
84 +  type(Properties), dimension(:),allocatable :: PropertyMap
85 +
86   contains
87  
88    subroutine setRlistDF( this_rlist )
# Line 55 | Line 93 | contains
93      rlistsq = rlist * rlist
94      
95      haveRlist = .true.
58 –    if( havePolicies ) do_forces_initialized = .true.
96  
97    end subroutine setRlistDF    
98  
99 +  subroutine createPropertyMap(status)
100 +    integer :: nAtypes
101 +    integer :: status
102 +    integer :: i
103 +    logical :: thisProperty
104 +    real (kind=DP) :: thisDPproperty
105 +
106 +    status = 0
107 +
108 +    nAtypes = getSize(atypes)
109 +
110 +    if (nAtypes == 0) then
111 +       status = -1
112 +       return
113 +    end if
114 +        
115 +    if (.not. allocated(PropertyMap)) then
116 +       allocate(PropertyMap(nAtypes))
117 +    endif
118 +
119 +    do i = 1, nAtypes
120 +       call getElementProperty(atypes, i, "is_LJ", thisProperty)
121 +       PropertyMap(i)%is_LJ = thisProperty
122 +
123 +       call getElementProperty(atypes, i, "is_Charge", thisProperty)
124 +       PropertyMap(i)%is_Charge = thisProperty
125 +      
126 +       if (thisProperty) then
127 +          call getElementProperty(atypes, i, "charge", thisDPproperty)
128 +          PropertyMap(i)%charge = thisDPproperty
129 +       endif
130 +
131 +       call getElementProperty(atypes, i, "is_DP", thisProperty)
132 +       PropertyMap(i)%is_DP = thisProperty
133 +
134 +       if (thisProperty) then
135 +          call getElementProperty(atypes, i, "dipole_moment", thisDPproperty)
136 +          PropertyMap(i)%dipole_moment = thisDPproperty
137 +       endif
138 +
139 +       call getElementProperty(atypes, i, "is_Sticky", thisProperty)
140 +       PropertyMap(i)%is_Sticky = thisProperty
141 +       call getElementProperty(atypes, i, "is_GB", thisProperty)
142 +       PropertyMap(i)%is_GB = thisProperty
143 +       call getElementProperty(atypes, i, "is_EAM", thisProperty)
144 +       PropertyMap(i)%is_EAM = thisProperty
145 +    end do
146 +
147 +    havePropertyMap = .true.
148 +
149 +  end subroutine createPropertyMap
150 +
151 +  subroutine setSimVariables()
152 +    SIM_uses_LJ = SimUsesLJ()
153 +    SIM_uses_sticky = SimUsesSticky()
154 +    SIM_uses_charges = SimUsesCharges()
155 +    SIM_uses_dipoles = SimUsesDipoles()
156 +    SIM_uses_RF = SimUsesRF()
157 +    SIM_uses_GB = SimUsesGB()
158 +    SIM_uses_EAM = SimUsesEAM()
159 +    SIM_requires_postpair_calc = SimRequiresPostpairCalc()
160 +    SIM_requires_prepair_calc = SimRequiresPrepairCalc()
161 +    SIM_uses_directional_atoms = SimUsesDirectionalAtoms()
162 +    SIM_uses_PBC = SimUsesPBC()
163 +
164 +    haveSIMvariables = .true.
165 +
166 +    return
167 +  end subroutine setSimVariables
168 +
169 +  subroutine doReadyCheck(error)
170 +    integer, intent(out) :: error
171 +
172 +    integer :: myStatus
173 +
174 +    error = 0
175 +    
176 +    if (.not. havePropertyMap) then
177 +
178 +       myStatus = 0
179 +
180 +       call createPropertyMap(myStatus)
181 +
182 +       if (myStatus .ne. 0) then
183 +          write(default_error, *) 'createPropertyMap failed in do_Forces!'
184 +          error = -1
185 +          return
186 +       endif
187 +    endif
188 +
189 +    if (.not. haveSIMvariables) then
190 +       call setSimVariables()
191 +    endif
192 +
193 +    if (.not. haveRlist) then
194 +       write(default_error, *) 'rList has not been set in do_Forces!'
195 +       error = -1
196 +       return
197 +    endif
198 +
199 +    if (SIM_uses_LJ .and. FF_uses_LJ) then
200 +       if (.not. havePolicies) then
201 +          write(default_error, *) 'LJ mixing Policies have not been set in do_Forces!'
202 +          error = -1
203 +          return
204 +       endif
205 +    endif
206 +
207 +    if (.not. haveNeighborList) then
208 +       write(default_error, *) 'neighbor list has not been initialized in do_Forces!'
209 +       error = -1
210 +       return
211 +    end if
212 +
213 +    if (.not. haveSaneForceField) then
214 +       write(default_error, *) 'Force Field is not sane in do_Forces!'
215 +       error = -1
216 +       return
217 +    end if
218 +
219 + #ifdef IS_MPI
220 +    if (.not. isMPISimSet()) then
221 +       write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
222 +       error = -1
223 +       return
224 +    endif
225 + #endif
226 +    return
227 +  end subroutine doReadyCheck
228 +    
229 +
230    subroutine init_FF(LJMIXPOLICY, use_RF_c, thisStat)
231  
232      integer, intent(in) :: LJMIXPOLICY
# Line 83 | Line 251 | contains
251    
252      FF_uses_LJ = .false.
253      FF_uses_sticky = .false.
254 +    FF_uses_charges = .false.
255      FF_uses_dipoles = .false.
256      FF_uses_GB = .false.
257      FF_uses_EAM = .false.
258      
259      call getMatchingElementList(atypes, "is_LJ", .true., nMatches, MatchList)
260      if (nMatches .gt. 0) FF_uses_LJ = .true.
261 +
262 +    call getMatchingElementList(atypes, "is_Charge", .true., nMatches, MatchList)
263 +    if (nMatches .gt. 0) FF_uses_charges = .true.
264      
265      call getMatchingElementList(atypes, "is_DP", .true., nMatches, MatchList)
266      if (nMatches .gt. 0) FF_uses_dipoles = .true.
# Line 103 | Line 275 | contains
275      call getMatchingElementList(atypes, "is_EAM", .true., nMatches, MatchList)
276      if (nMatches .gt. 0) FF_uses_EAM = .true.
277      
278 +    !! Assume sanity (for the sake of argument)
279 +    haveSaneForceField = .true.
280 +
281      !! check to make sure the FF_uses_RF setting makes sense
282      
283      if (FF_uses_dipoles) then
# Line 114 | Line 289 | contains
289         if (FF_uses_RF) then          
290            write(default_error,*) 'Using Reaction Field with no dipoles?  Huh?'
291            thisStat = -1
292 +          haveSaneForceField = .false.
293            return
294         endif
295      endif
# Line 128 | Line 304 | contains
304         case default
305            write(default_error,*) 'unknown LJ Mixing Policy!'
306            thisStat = -1
307 +          haveSaneForceField = .false.
308            return            
309         end select
310         if (my_status /= 0) then
311            thisStat = -1
312 +          haveSaneForceField = .false.
313            return
314         end if
315 +       havePolicies = .true.
316      endif
317  
318      if (FF_uses_sticky) then
319         call check_sticky_FF(my_status)
320         if (my_status /= 0) then
321            thisStat = -1
322 +          haveSaneForceField = .false.
323            return
324         end if
325      endif
326  
327  
328      if (FF_uses_EAM) then
329 <       call init_EAM_FF(my_status)
329 >         call init_EAM_FF(my_status)
330         if (my_status /= 0) then
331 +          write(default_error, *) "init_EAM_FF returned a bad status"
332            thisStat = -1
333 +          haveSaneForceField = .false.
334            return
335         end if
336      endif
337  
156 –
157 –    
338      if (FF_uses_GB) then
339         call check_gb_pair_FF(my_status)
340         if (my_status .ne. 0) then
341            thisStat = -1
342 +          haveSaneForceField = .false.
343            return
344         endif
345      endif
346  
347      if (FF_uses_GB .and. FF_uses_LJ) then
348      endif
349 <    if (.not. do_forces_initialized) then
349 >    if (.not. haveNeighborList) then
350         !! Create neighbor lists
351 <       call expandNeighborList(getNlocal(), my_status)
351 >       call expandNeighborList(nLocal, my_status)
352         if (my_Status /= 0) then
353            write(default_error,*) "SimSetup: ExpandNeighborList returned error."
354            thisStat = -1
355            return
356         endif
357 +       haveNeighborList = .true.
358      endif
359      
178 –
179 –    havePolicies = .true.
180 –    if( haveRlist ) do_forces_initialized = .true.
181 –
360    end subroutine init_FF
361    
362  
# Line 187 | Line 365 | contains
365    subroutine do_force_loop(q, A, u_l, f, t, tau, pot, do_pot_c, do_stress_c, &
366         error)
367      !! Position array provided by C, dimensioned by getNlocal
368 <    real ( kind = dp ), dimension(3,getNlocal()) :: q
368 >    real ( kind = dp ), dimension(3,nLocal) :: q
369      !! Rotation Matrix for each long range particle in simulation.
370 <    real( kind = dp), dimension(9,getNlocal()) :: A    
370 >    real( kind = dp), dimension(9,nLocal) :: A    
371      !! Unit vectors for dipoles (lab frame)
372 <    real( kind = dp ), dimension(3,getNlocal()) :: u_l
372 >    real( kind = dp ), dimension(3,nLocal) :: u_l
373      !! Force array provided by C, dimensioned by getNlocal
374 <    real ( kind = dp ), dimension(3,getNlocal()) :: f
374 >    real ( kind = dp ), dimension(3,nLocal) :: f
375      !! Torsion array provided by C, dimensioned by getNlocal
376 <    real( kind = dp ), dimension(3,getNlocal()) :: t    
376 >    real( kind = dp ), dimension(3,nLocal) :: t    
377 >
378      !! Stress Tensor
379      real( kind = dp), dimension(9) :: tau  
380      real ( kind = dp ) :: pot
# Line 206 | Line 385 | contains
385      real( kind = DP ) :: pot_local
386      integer :: nrow
387      integer :: ncol
388 +    integer :: nprocs
389   #endif
210 –    integer :: nlocal
390      integer :: natoms    
391      logical :: update_nlist  
392      integer :: i, j, jbeg, jend, jnab
# Line 215 | Line 394 | contains
394      real( kind = DP ) ::  rijsq
395      real(kind=dp),dimension(3) :: d
396      real(kind=dp) :: rfpot, mu_i, virial
397 <    integer :: me_i
397 >    integer :: me_i, me_j
398      logical :: is_dp_i
399      integer :: neighborListSize
400      integer :: listerror, error
401      integer :: localError
402 +    integer :: propPack_i, propPack_j
403  
404 <    real(kind=dp) :: listSkin = 1.0
225 <    
404 >    real(kind=dp) :: listSkin = 1.0  
405  
406      !! initialize local variables  
407  
408   #ifdef IS_MPI
409      pot_local = 0.0_dp
231 –    nlocal = getNlocal()
410      nrow   = getNrow(plan_row)
411      ncol   = getNcol(plan_col)
412   #else
235 –    nlocal = getNlocal()
413      natoms = nlocal
414   #endif
415  
416 <    call check_initialization(localError)
416 >    call doReadyCheck(localError)
417      if ( localError .ne. 0 ) then
418 <       call handleError("do_force_loop","Not Initialized")
418 >       call handleError("do_force_loop", "Not Initialized")
419         error = -1
420         return
421      end if
# Line 247 | Line 424 | contains
424      do_pot = do_pot_c
425      do_stress = do_stress_c
426  
250 –
427      ! Gather all information needed by all force loops:
428      
429   #ifdef IS_MPI    
# Line 255 | Line 431 | contains
431      call gather(q,q_Row,plan_row3d)
432      call gather(q,q_Col,plan_col3d)
433          
434 <    if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
434 >    if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
435         call gather(u_l,u_l_Row,plan_row3d)
436         call gather(u_l,u_l_Col,plan_col3d)
437        
# Line 264 | Line 440 | contains
440      endif
441      
442   #endif
443 +
444 + !! Begin force loop timing:
445 + #ifdef PROFILE
446 +    call cpu_time(forceTimeInitial)
447 +    nloops = nloops + 1
448 + #endif
449    
450 <    if (FF_RequiresPrepairCalc() .and. SimRequiresPrepairCalc()) then
450 >    if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
451         !! See if we need to update neighbor lists
452         call checkNeighborList(nlocal, q, listSkin, update_nlist)  
453         !! if_mpi_gather_stuff_for_prepair
454         !! do_prepair_loop_if_needed
455         !! if_mpi_scatter_stuff_from_prepair
456         !! if_mpi_gather_stuff_from_prepair_to_main_loop
457 <
457 >    
458   !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
459   #ifdef IS_MPI
460      
# Line 348 | Line 530 | contains
530         neighborListSize = size(list)
531    
532         nlist = 0
533 <      
533 >
534         do i = 1, natoms-1
535            point(i) = nlist + 1
536            
# Line 360 | Line 542 | contains
542            
543  
544               if (rijsq < rlistsq) then
545 <                
545 >
546 >          
547                  nlist = nlist + 1
548                
549                  if (nlist > neighborListSize) then
# Line 385 | Line 568 | contains
568         point(natoms) = nlist + 1
569        
570      else !! (update)
571 <      
571 >  
572         ! use the list to find the neighbors
573         do i = 1, natoms-1
574            JBEG = POINT(i)
# Line 406 | Line 589 | contains
589      endif    
590   #endif
591      !! Do rest of preforce calculations
592 <   call do_preforce(nlocal,pot)
592 >    !! do necessary preforce calculations  
593 >    call do_preforce(nlocal,pot)
594 >   ! we have already updated the neighbor list set it to false...
595 >   update_nlist = .false.
596      else
597         !! See if we need to update neighbor lists for non pre-pair
598         call checkNeighborList(nlocal, q, listSkin, update_nlist)  
# Line 425 | Line 611 | contains
611   #ifdef IS_MPI
612      
613      if (update_nlist) then
428 –      
614         !! save current configuration, construct neighbor list,
615         !! and calculate forces
616         call saveNeighborList(nlocal, q)
# Line 434 | Line 619 | contains
619         nlist = 0      
620        
621         do i = 1, nrow
622 +
623            point(i) = nlist + 1
624            
625            inner: do j = 1, ncol
# Line 491 | Line 677 | contains
677   #else
678      
679      if (update_nlist) then
680 <      
680 >
681         ! save current configuration, contruct neighbor list,
682         ! and calculate forces
683         call saveNeighborList(natoms, q)
# Line 559 | Line 745 | contains
745   #endif
746      
747      ! phew, done with main loop.
748 <    
748 >
749 > !! Do timing
750 > #ifdef PROFILE
751 >    call cpu_time(forceTimeFinal)
752 >    forceTime = forceTime + forceTimeFinal - forceTimeInitial
753 > #endif
754 >
755 >
756   #ifdef IS_MPI
757      !!distribute forces
758    
# Line 575 | Line 768 | contains
768         f(1:3,i) = f(1:3,i) + f_temp(1:3,i)
769      end do
770      
771 <    if (FF_UsesDirectionalAtoms() .and. SimUsesDirectionalAtoms()) then
771 >    if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
772         t_temp = 0.0_dp
773         call scatter(t_Row,t_temp,plan_row3d)
774         do i = 1,nlocal
# Line 609 | Line 802 | contains
802      endif    
803   #endif
804  
805 <    if (FF_RequiresPostpairCalc() .and. SimRequiresPostpairCalc()) then
805 >    if (FF_RequiresPostpairCalc() .and. SIM_requires_postpair_calc) then
806        
807 <       if (FF_uses_RF .and. SimUsesRF()) then
807 >       if (FF_uses_RF .and. SIM_uses_RF) then
808            
809   #ifdef IS_MPI
810            call scatter(rf_Row,rf,plan_row3d)
# Line 621 | Line 814 | contains
814            end do
815   #endif
816            
817 <          do i = 1, getNlocal()
817 >          do i = 1, nLocal
818  
819               rfpot = 0.0_DP
820   #ifdef IS_MPI
# Line 629 | Line 822 | contains
822   #else
823               me_i = atid(i)
824   #endif
825 <             call getElementProperty(atypes, me_i, "is_DP", is_DP_i)      
826 <             if ( is_DP_i ) then
827 <                call getElementProperty(atypes, me_i, "dipole_moment", mu_i)
825 >
826 >             if (PropertyMap(me_i)%is_DP) then
827 >
828 >                mu_i = PropertyMap(me_i)%dipole_moment
829 >
830                  !! The reaction field needs to include a self contribution
831                  !! to the field:
832 <                call accumulate_self_rf(i, mu_i, u_l)            
832 >                call accumulate_self_rf(i, mu_i, u_l)
833                  !! Get the reaction field contribution to the
834                  !! potential and torques:
835                  call reaction_field_final(i, mu_i, u_l, rfpot, t, do_pot)
# Line 671 | Line 866 | contains
866         tau = tau_Temp
867         virial = virial_Temp
868      endif
869 <
869 >    
870   #endif
871      
872 +    
873 +    
874    end subroutine do_force_loop
875  
876    subroutine do_pair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
877  
878      real( kind = dp ) :: pot
879 <    real( kind = dp ), dimension(3,getNlocal()) :: u_l
880 <    real (kind=dp), dimension(9,getNlocal()) :: A
881 <    real (kind=dp), dimension(3,getNlocal()) :: f
882 <    real (kind=dp), dimension(3,getNlocal()) :: t
879 >    real( kind = dp ), dimension(3,nLocal) :: u_l
880 >    real (kind=dp), dimension(9,nLocal) :: A
881 >    real (kind=dp), dimension(3,nLocal) :: f
882 >    real (kind=dp), dimension(3,nLocal) :: t
883  
884      logical, intent(inout) :: do_pot, do_stress
885      integer, intent(in) :: i, j
# Line 695 | Line 892 | contains
892      logical :: is_EAM_i,is_EAM_j
893      logical :: is_Sticky_i, is_Sticky_j
894      integer :: me_i, me_j
895 <
895 >    integer :: propPack_i
896 >    integer :: propPack_j
897      r = sqrt(rijsq)
898  
899   #ifdef IS_MPI
# Line 712 | Line 910 | contains
910      me_j = atid(j)
911  
912   #endif
913 +    
914 +    if (FF_uses_LJ .and. SIM_uses_LJ) then
915  
916 <    if (FF_uses_LJ .and. SimUsesLJ()) then
717 <       call getElementProperty(atypes, me_i, "is_LJ", is_LJ_i)
718 <       call getElementProperty(atypes, me_j, "is_LJ", is_LJ_j)
719 <
720 <       if ( is_LJ_i .and. is_LJ_j ) &
916 >       if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) &
917              call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
918 +
919      endif
920  
921 <    if (FF_uses_dipoles .and. SimUsesDipoles()) then
725 <       call getElementProperty(atypes, me_i, "is_DP", is_DP_i)
726 <       call getElementProperty(atypes, me_j, "is_DP", is_DP_j)
921 >    if (FF_uses_dipoles .and. SIM_uses_dipoles) then
922        
923 <       if ( is_DP_i .and. is_DP_j ) then
729 <          
923 >       if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
924            call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
925                 do_pot, do_stress)
926 <          if (FF_uses_RF .and. SimUsesRF()) then
926 >          if (FF_uses_RF .and. SIM_uses_RF) then
927               call accumulate_rf(i, j, r, u_l)
928               call rf_correct_forces(i, j, d, r, u_l, f, do_stress)
929            endif
# Line 737 | Line 931 | contains
931         endif
932      endif
933  
934 <    if (FF_uses_Sticky .and. SimUsesSticky()) then
934 >    if (FF_uses_Sticky .and. SIM_uses_sticky) then
935  
936 <       call getElementProperty(atypes, me_i, "is_Sticky", is_Sticky_i)
743 <       call getElementProperty(atypes, me_j, "is_Sticky", is_Sticky_j)
744 <
745 <       if ( is_Sticky_i .and. is_Sticky_j ) then
936 >       if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
937            call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
938                 do_pot, do_stress)
939         endif
940      endif
941  
942  
943 <    if (FF_uses_GB .and. SimUsesGB()) then
753 <
754 <       call getElementProperty(atypes, me_i, "is_GB", is_GB_i)
755 <       call getElementProperty(atypes, me_j, "is_GB", is_GB_j)
943 >    if (FF_uses_GB .and. SIM_uses_GB) then
944        
945 <       if ( is_GB_i .and. is_GB_j ) then
945 >       if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
946            call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
947                 do_pot, do_stress)          
948         endif
949 +
950      endif
951      
952  
953    
954 <   if (FF_uses_EAM .and. SimUsesEAM()) then
766 <      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
767 <      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
954 >   if (FF_uses_EAM .and. SIM_uses_EAM) then
955        
956 <      if ( is_EAM_i .and. is_EAM_j ) &
957 <           call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
956 >      if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
957 >         call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
958 >      endif
959 >
960     endif
772 –
961  
774 –
775 –
962    end subroutine do_pair
963  
964  
965  
966    subroutine do_prepair(i, j, rijsq, d, do_pot, do_stress, u_l, A, f, t, pot)
967     real( kind = dp ) :: pot
968 <   real( kind = dp ), dimension(3,getNlocal()) :: u_l
969 <   real (kind=dp), dimension(9,getNlocal()) :: A
970 <   real (kind=dp), dimension(3,getNlocal()) :: f
971 <   real (kind=dp), dimension(3,getNlocal()) :: t
968 >   real( kind = dp ), dimension(3,nLocal) :: u_l
969 >   real (kind=dp), dimension(9,nLocal) :: A
970 >   real (kind=dp), dimension(3,nLocal) :: f
971 >   real (kind=dp), dimension(3,nLocal) :: t
972    
973     logical, intent(inout) :: do_pot, do_stress
974     integer, intent(in) :: i, j
# Line 811 | Line 997 | contains
997     me_j = atid(j)
998    
999   #endif
1000 <  
1001 <   if (FF_uses_EAM .and. SimUsesEAM()) then
1002 <      call getElementProperty(atypes, me_i, "is_EAM", is_EAM_i)
1003 <      call getElementProperty(atypes, me_j, "is_EAM", is_EAM_j)
818 <      
819 <      if ( is_EAM_i .and. is_EAM_j ) &
1000 >    
1001 >   if (FF_uses_EAM .and. SIM_uses_EAM) then
1002 >
1003 >      if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
1004             call calc_EAM_prepair_rho(i, j, d, r, rijsq )
1005 +
1006     endif
1007 <  end subroutine do_prepair
1007 >  
1008 > end subroutine do_prepair
1009  
1010  
1011  
# Line 828 | Line 1014 | contains
1014      integer :: nlocal
1015      real( kind = dp ) :: pot
1016  
1017 <    if (FF_uses_EAM .and. SimUsesEAM()) then
1017 >    if (FF_uses_EAM .and. SIM_uses_EAM) then
1018         call calc_EAM_preforce_Frho(nlocal,pot)
1019      endif
1020  
# Line 847 | Line 1033 | contains
1033      d(1:3) = q_j(1:3) - q_i(1:3)
1034  
1035      ! Wrap back into periodic box if necessary
1036 <    if ( SimUsesPBC() ) then
1036 >    if ( SIM_uses_PBC ) then
1037        
1038         if( .not.boxIsOrthorhombic ) then
1039            ! calc the scaled coordinates.
# Line 887 | Line 1073 | contains
1073      
1074    end subroutine get_interatomic_vector
1075    
890 –  subroutine check_initialization(error)
891 –    integer, intent(out) :: error
892 –    
893 –    error = 0
894 –    ! Make sure we are properly initialized.
895 –    if (.not. do_forces_initialized) then
896 –       write(*,*) "Forces not initialized"
897 –       error = -1
898 –       return
899 –    endif
900 –
901 – #ifdef IS_MPI
902 –    if (.not. isMPISimSet()) then
903 –       write(default_error,*) "ERROR: mpiSimulation has not been initialized!"
904 –       error = -1
905 –       return
906 –    endif
907 – #endif
908 –    
909 –    return
910 –  end subroutine check_initialization
911 –
912 –  
1076    subroutine zero_work_arrays()
1077      
1078   #ifdef IS_MPI
# Line 941 | Line 1104 | contains
1104  
1105   #endif
1106  
1107 +
1108 +    if (FF_uses_EAM .and. SIM_uses_EAM) then
1109 +       call clean_EAM()
1110 +    endif
1111 +
1112 +
1113 +
1114 +
1115 +
1116      rf = 0.0_dp
1117      tau_Temp = 0.0_dp
1118      virial_Temp = 0.0_dp
# Line 1053 | Line 1225 | end module do_Forces
1225      doesit = FF_uses_RF
1226    end function FF_RequiresPostpairCalc
1227    
1228 + #ifdef PROFILE
1229 +  function getforcetime() result(totalforcetime)
1230 +    real(kind=dp) :: totalforcetime
1231 +    totalforcetime = forcetime
1232 +  end function getforcetime
1233 + #endif
1234 +
1235 + !! This cleans componets of force arrays belonging only to fortran
1236 +
1237   end module do_Forces

Diff Legend

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