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 1160 by gezelter, Tue May 11 21:31:15 2004 UTC vs.
Revision 1197 by gezelter, Wed May 26 16:41:23 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.55 2004-05-11 21:31:14 gezelter Exp $, $Date: 2004-05-11 21:31:14 $, $Name: not supported by cvs2svn $, $Revision: 1.55 $
7 > !! @version $Id: do_Forces.F90,v 1.62 2004-05-26 16:41:23 gezelter Exp $, $Date: 2004-05-26 16:41:23 $, $Name: not supported by cvs2svn $, $Revision: 1.62 $
8  
9   module do_Forces
10    use force_globals
# Line 32 | Line 32 | module do_Forces
32   #define __FORTRAN90
33   #include "fForceField.h"
34   #include "fSwitchingFunction.h"
35 +
36 +  INTEGER, PARAMETER:: PREPAIR_LOOP = 1
37 +  INTEGER, PARAMETER:: PAIR_LOOP    = 2
38  
39    logical, save :: haveRlist = .false.
40    logical, save :: haveNeighborList = .false.
# Line 370 | Line 373 | contains
373    subroutine do_force_loop(q, q_group, A, u_l, f, t, tau, pot, &
374         do_pot_c, do_stress_c, error)
375      !! Position array provided by C, dimensioned by getNlocal
376 <    real ( kind = dp ), dimension(3,nLocal) :: q
377 <    !! molecular center-of-mass position array
378 <    real ( kind = dp ), dimension(3,nGroup) :: q_group
376 >    real ( kind = dp ), dimension(3, nLocal) :: q
377 >    !! molecular center-of-mass position array
378 >    real ( kind = dp ), dimension(3, nGroup) :: q_group
379      !! Rotation Matrix for each long range particle in simulation.
380 <    real( kind = dp), dimension(9,nLocal) :: A    
380 >    real( kind = dp), dimension(9, nLocal) :: A    
381      !! Unit vectors for dipoles (lab frame)
382      real( kind = dp ), dimension(3,nLocal) :: u_l
383      !! Force array provided by C, dimensioned by getNlocal
# Line 399 | Line 402 | contains
402   #endif
403      integer :: natoms    
404      logical :: update_nlist  
405 <    integer :: i, j, jbeg, jend, jnab
405 >    integer :: i, j, jstart, jend, jnab
406 >    integer :: istart, iend
407      integer :: ia, jb, atom1, atom2
408      integer :: nlist
409 <    real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vab
409 >    real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vij
410      real( kind = DP ) :: sw, dswdr, swderiv, mf
411 <    real(kind=dp),dimension(3) :: d_atm, d_grp
411 >    real(kind=dp),dimension(3) :: d_atm, d_grp, fpair, fij
412      real(kind=dp) :: rfpot, mu_i, virial
413 <    integer :: me_i, me_j
413 >    integer :: me_i, me_j, n_in_i, n_in_j
414      logical :: is_dp_i
415      integer :: neighborListSize
416      integer :: listerror, error
417      integer :: localError
418      integer :: propPack_i, propPack_j
419 +    integer :: loopStart, loopEnd, loop
420  
421      real(kind=dp) :: listSkin = 1.0  
422      
# Line 464 | Line 469 | contains
469      nloops = nloops + 1
470   #endif
471      
472 +    loopEnd = PAIR_LOOP
473      if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
474 <       !! See if we need to update neighbor lists
475 <
476 <       call checkNeighborList(nGroup, q_group, listSkin, update_nlist)
474 >       loopStart = PREPAIR_LOOP
475 >    else
476 >       loopStart = PAIR_LOOP
477 >    endif
478  
479 <       !! if_mpi_gather_stuff_for_prepair
480 <       !! do_prepair_loop_if_needed
481 <       !! if_mpi_scatter_stuff_from_prepair
482 <       !! if_mpi_gather_stuff_from_prepair_to_main_loop
479 >    do loop = loopStart, loopEnd
480 >
481 >       ! See if we need to update neighbor lists
482 >       ! (but only on the first time through):
483 >       if (loop .eq. loopStart) then
484 >          call checkNeighborList(nGroup, q_group, listSkin, update_nlist)
485 >       endif
486        
477 –       !--------------------PREFORCE LOOP----------->>>>>>>>>>>>>>>>>>>>>>>>>>>
478 – #ifdef IS_MPI
479 –      
487         if (update_nlist) then
488 <          
489 <          !! save current configuration, construct neighbor list,
483 <          !! and calculate forces
484 <
485 <          call saveNeighborList(nGroup, q_group)
486 <          
488 >          !! save current configuration and construct neighbor list
489 >          call saveNeighborList(nGroup, q_group)          
490            neighborListSize = size(list)
491 <          nlist = 0      
489 <          
490 <          do i = 1, nrow_group
491 <             point(i) = nlist + 1
492 <            
493 <             do j = 1, ncol_group
494 <                
495 <                call get_interatomic_vector(q_group_Row(:,i), &
496 <                     q_group_Col(:,j), d_grp, rgrpsq)
497 <                
498 <                if (rgrpsq < rlistsq) then
499 <                   nlist = nlist + 1
500 <                  
501 <                   if (nlist > neighborListSize) then
502 <                      call expandNeighborList(nGroup, listerror)
503 <                      if (listerror /= 0) then
504 <                         error = -1
505 <                         write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
506 <                         return
507 <                      end if
508 <                      neighborListSize = size(list)
509 <                   endif
510 <                  
511 <                   list(nlist) = j
512 <                  
513 <                   do ia = groupStart(i), groupStart(i+1)-1
514 <                      atom1 = groupList(ia)
515 <
516 <                      prepair_inner: do jb = groupStart(j), groupStart(j+1)-1
517 <                         atom2 = groupList(jb)
518 <                        
519 <                         if (skipThisPair(atom1, atom2)) cycle prepair_inner
520 <                        
521 <                         call get_interatomic_vector(q_Row(:,atom1), &
522 <                              q_Col(:,atom2), d_atm, ratmsq)
523 <                        
524 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, &
525 <                              rgrpsq, d_grp, do_pot, do_stress, &
526 <                              u_l, A, f, t, pot_local)
527 <
528 <                      enddo prepair_inner
529 <                   enddo
530 <                end if
531 <             enddo
532 <          enddo
533 <          point(nrow_group + 1) = nlist + 1          
534 <                          
535 <       else  !! (of update_check)
536 <          
537 <          ! use the list to find the neighbors
538 <          do i = 1, nrow_group
539 <             JBEG = POINT(i)
540 <             JEND = POINT(i+1) - 1
541 <             ! check that group i has neighbors
542 <             if (jbeg .le. jend) then
543 <                
544 <                do jnab = jbeg, jend
545 <                   j = list(jnab)
546 <                  
547 <                   do ia = groupStart(i), groupStart(i+1)-1
548 <                      atom1 = groupList(ia)
549 <
550 <                      do jb = groupStart(j), groupStart(j+1)-1
551 <                         atom2 = groupList(jb)                        
552 <                        
553 <                         call get_interatomic_vector(q_Row(:,atom1), &
554 <                              q_Col(:,atom2), d_atm, ratmsq)
555 <                        
556 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, &
557 <                              rgrpsq, d_grp, do_pot, do_stress, &
558 <                              u_l, A, f, t, pot_local)
559 <                        
560 <                      enddo
561 <                   enddo
562 <                enddo
563 <             endif
564 <          enddo
491 >          nlist = 0
492         endif
493        
494 +       istart = 1
495 + #ifdef IS_MPI
496 +       iend = nrow_group
497   #else
498 <
499 <       if (update_nlist) then
498 >       iend = nGroup - 1
499 > #endif
500 >       outer: do i = istart, iend
501            
502 <          !! save current configuration, construct neighbor list,
572 <          !! and calculate forces
502 >          if (update_nlist) point(i) = nlist + 1
503            
504 <          call saveNeighborList(nGroup, q_group)
504 >          n_in_i = groupStart(i+1) - groupStart(i)
505            
506 <          neighborListSize = size(list)
507 <          nlist = 0      
506 >          if (update_nlist) then
507 > #ifdef IS_MPI
508 >             jstart = 1
509 >             jend = ncol_group
510 > #else
511 >             jstart = i+1
512 >             jend = nGroup
513 > #endif
514 >          else            
515 >             jstart = point(i)
516 >             jend = point(i+1) - 1
517 >             ! make sure group i has neighbors
518 >             if (jstart .gt. jend) cycle outer
519 >          endif
520            
521 <          do i = 1, nGroup-1
522 <             point(i) = nlist + 1
523 <            
524 <             do j = i+1, nGroup
525 <                
526 <                call get_interatomic_vector(q_group(:,i), &
527 <                     q_group(:,j), d_grp, rgrpsq)
528 <                
529 <                if (rgrpsq < rlistsq) then
521 >          do jnab = jstart, jend
522 >             if (update_nlist) then
523 >                j = jnab
524 >             else
525 >                j = list(jnab)
526 >             endif
527 > #ifdef IS_MPI
528 >             call get_interatomic_vector(q_group_Row(:,i), &
529 >                  q_group_Col(:,j), d_grp, rgrpsq)
530 > #else
531 >             call get_interatomic_vector(q_group(:,i), &
532 >                  q_group(:,j), d_grp, rgrpsq)
533 > #endif
534 >             if (rgrpsq < rlistsq) then
535 >                if (update_nlist) then
536                     nlist = nlist + 1
537                    
538                     if (nlist > neighborListSize) then
# Line 598 | Line 546 | contains
546                     endif
547                    
548                     list(nlist) = j
549 <                  
602 <                   do ia = groupStart(i), groupStart(i+1)-1
603 <                      atom1 = groupList(ia)
604 <                      
605 <                      prepair_inner: do jb = groupStart(j), groupStart(j+1)-1
606 <                         atom2 = groupList(jb)
607 <                        
608 <                         if (skipThisPair(atom1, atom2)) cycle prepair_inner
609 <                        
610 <                         call get_interatomic_vector(q(:,atom1), &
611 <                              q(:,atom2), d_atm, ratmsq)
612 <                        
613 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, &
614 <                              rgrpsq, d_grp, do_pot, do_stress, &
615 <                              u_l, A, f, t, pot)
616 <                        
617 <                      enddo prepair_inner
618 <                   enddo
619 <                end if
620 <             enddo
621 <          enddo
622 <          point(nGroup) = nlist + 1          
623 <                          
624 <       else  !! (of update_check)
625 <          
626 <          ! use the list to find the neighbors
627 <          do i = 1, nGroup-1
628 <             JBEG = POINT(i)
629 <             JEND = POINT(i+1) - 1
630 <             ! check that group i has neighbors
631 <             if (jbeg .le. jend) then
549 >                endif
550                  
551 <                do jnab = jbeg, jend
552 <                   j = list(jnab)
553 <                  
636 <                   do ia = groupStart(i), groupStart(i+1)-1
637 <                      atom1 = groupList(ia)
638 <
639 <                      do jb = groupStart(j), groupStart(j+1)-1
640 <                         atom2 = groupList(jb)                        
641 <                        
642 <                         call get_interatomic_vector(q(:,atom1), &
643 <                              q(:,atom2), d_atm, ratmsq)
644 <                        
645 <                         call do_prepair(atom1, atom2, ratmsq, d_atm, &
646 <                              rgrpsq, d_grp, do_pot, do_stress, &
647 <                              u_l, A, f, t, pot)
648 <                        
649 <                      enddo
650 <                   enddo
651 <                enddo
652 <             endif
653 <          enddo
654 <       endif
655 <      
656 < #endif
657 <      
658 <       !! Do rest of preforce calculations
659 <       !! do necessary preforce calculations  
660 <       call do_preforce(nlocal,pot)
661 <       ! we have already updated the neighbor list set it to false...
662 <       update_nlist = .false.
663 <    else
664 <       !! See if we need to update neighbor lists for non pre-pair
665 <       call checkNeighborList(nGroup, q_group, listSkin, update_nlist)  
666 <    endif
667 <    
668 <    !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>
669 <    
670 < #ifdef IS_MPI
671 <    
672 <    if (update_nlist) then
673 <      
674 <       !! save current configuration, construct neighbor list,
675 <       !! and calculate forces
676 <      
677 <       call saveNeighborList(nGroup, q_group)
678 <      
679 <       neighborListSize = size(list)
680 <       nlist = 0      
681 <      
682 <       do i = 1, nrow_group
683 <          point(i) = nlist + 1
684 <          
685 <          do j = 1, ncol_group
686 <            
687 <             call get_interatomic_vector(q_group_Row(:,i), &
688 <                  q_group_Col(:,j), d_grp, rgrpsq)
689 <            
690 <             if (rgrpsq < rlistsq) then
691 <                nlist = nlist + 1
692 <                
693 <                if (nlist > neighborListSize) then
694 <                   call expandNeighborList(nGroup, listerror)
695 <                   if (listerror /= 0) then
696 <                      error = -1
697 <                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
698 <                      return
699 <                   end if
700 <                   neighborListSize = size(list)
551 >                if (loop .eq. PAIR_LOOP) then
552 >                   vij = 0.0d0
553 >                   fij(1:3) = 0.0d0
554                  endif
555                  
703 –                list(nlist) = j
704 –                
705 –                vab = 0.0d0
556                  call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
557                       in_switching_region)
558                  
559 +                n_in_j = groupStart(j+1) - groupStart(j)
560 +                
561                  do ia = groupStart(i), groupStart(i+1)-1
562 +                  
563                     atom1 = groupList(ia)
564                    
565                     inner: do jb = groupStart(j), groupStart(j+1)-1
566 +                      
567                        atom2 = groupList(jb)
568                        
569                        if (skipThisPair(atom1, atom2)) cycle inner
570                        
571 <                      call get_interatomic_vector(q_Row(:,atom1), &
572 <                           q_Col(:,atom2), d_atm, ratmsq)
573 <                      
574 <                      call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
575 <                           do_pot, do_stress, &
576 <                           u_l, A, f, t, pot_local, vpair)
577 <                      
578 <                      vab = vab + vpair
571 >                      if ((n_in_i .eq. 1).and.(n_in_j .eq. 1)) then
572 >                         d_atm(1:3) = d_grp(1:3)
573 >                         ratmsq = rgrpsq
574 >                      else
575 > #ifdef IS_MPI
576 >                         call get_interatomic_vector(q_Row(:,atom1), &
577 >                              q_Col(:,atom2), d_atm, ratmsq)
578 > #else
579 >                         call get_interatomic_vector(q(:,atom1), &
580 >                              q(:,atom2), d_atm, ratmsq)
581 > #endif
582 >                      endif
583 >                      if (loop .eq. PREPAIR_LOOP) then
584 > #ifdef IS_MPI                      
585 >                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
586 >                              rgrpsq, d_grp, do_pot, do_stress, &
587 >                              u_l, A, f, t, pot_local)
588 > #else
589 >                         call do_prepair(atom1, atom2, ratmsq, d_atm, sw, &
590 >                              rgrpsq, d_grp, do_pot, do_stress, &
591 >                              u_l, A, f, t, pot)
592 > #endif                                              
593 >                      else
594 > #ifdef IS_MPI                      
595 >                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
596 >                              do_pot, &
597 >                              u_l, A, f, t, pot_local, vpair, fpair)
598 > #else
599 >                         call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
600 >                              do_pot,  &
601 >                              u_l, A, f, t, pot, vpair, fpair)
602 > #endif
603 >                         vij = vij + vpair
604 >                         fij(1:3) = fij(1:3) + fpair(1:3)
605 >                      endif
606                     enddo inner
607                  enddo
608 <
609 <                if (in_switching_region) then
610 <                   swderiv = vab*dswdr/rgrp
611 <                  
612 <                   do ia=groupStart(i), groupStart(i+1)-1
613 <                      atom1=groupList(ia)
614 <                      mf = mfact(atom1)                  
615 <                      f_Row(1,atom1) = f_Row(1,atom1) - swderiv*d_grp(1)*mf
616 <                      f_Row(2,atom1) = f_Row(2,atom1) - swderiv*d_grp(2)*mf
617 <                      f_Row(3,atom1) = f_Row(3,atom1) - swderiv*d_grp(3)*mf
618 <                   enddo
608 >                
609 >                if (loop .eq. PAIR_LOOP) then
610 >                   if (in_switching_region) then
611 >                      swderiv = vij*dswdr/rgrp
612 >                      fij(1) = fij(1) + swderiv*d_grp(1)
613 >                      fij(2) = fij(2) + swderiv*d_grp(2)
614 >                      fij(3) = fij(3) + swderiv*d_grp(3)
615 >                      
616 >                      do ia=groupStart(i), groupStart(i+1)-1
617 >                         atom1=groupList(ia)
618 >                         mf = mfact(atom1)
619 > #ifdef IS_MPI
620 >                         f_Row(1,atom1) = f_Row(1,atom1) + swderiv*d_grp(1)*mf
621 >                         f_Row(2,atom1) = f_Row(2,atom1) + swderiv*d_grp(2)*mf
622 >                         f_Row(3,atom1) = f_Row(3,atom1) + swderiv*d_grp(3)*mf
623 > #else
624 >                         f(1,atom1) = f(1,atom1) + swderiv*d_grp(1)*mf
625 >                         f(2,atom1) = f(2,atom1) + swderiv*d_grp(2)*mf
626 >                         f(3,atom1) = f(3,atom1) + swderiv*d_grp(3)*mf
627 > #endif
628 >                      enddo
629 >                      
630 >                      do jb=groupStart(j), groupStart(j+1)-1
631 >                         atom2=groupList(jb)
632 >                         mf = mfact(atom2)
633 > #ifdef IS_MPI
634 >                         f_Col(1,atom2) = f_Col(1,atom2) - swderiv*d_grp(1)*mf
635 >                         f_Col(2,atom2) = f_Col(2,atom2) - swderiv*d_grp(2)*mf
636 >                         f_Col(3,atom2) = f_Col(3,atom2) - swderiv*d_grp(3)*mf
637 > #else
638 >                         f(1,atom2) = f(1,atom2) - swderiv*d_grp(1)*mf
639 >                         f(2,atom2) = f(2,atom2) - swderiv*d_grp(2)*mf
640 >                         f(3,atom2) = f(3,atom2) - swderiv*d_grp(3)*mf
641 > #endif
642 >                      enddo
643 >                   endif
644                    
645 <                   do jb=groupStart(j), groupStart(j+1)-1
740 <                      atom2=groupList(jb)
741 <                      mf = mfact(atom2)
742 <                      f_Col(1,atom2) = f_Col(1,atom2) + swderiv*d_grp(1)*mf
743 <                      f_Col(2,atom2) = f_Col(2,atom2) + swderiv*d_grp(2)*mf
744 <                      f_Col(3,atom2) = f_Col(3,atom2) + swderiv*d_grp(3)*mf
745 <                   enddo                        
645 >                   if (do_stress) call add_stress_tensor(d_grp, fij)
646                  endif
747 –
647               end if
648            enddo
649 <       enddo
751 <       point(nrow_group + 1) = nlist + 1          
649 >       enddo outer
650        
651 <    else  !! (of update_check)
652 <      
653 <       ! use the list to find the neighbors
654 <       do i = 1, nrow_group
655 <          JBEG = POINT(i)
656 <          JEND = POINT(i+1) - 1
657 <          ! check that group i has neighbors
658 <          if (jbeg .le. jend) then
659 <            
660 <             do jnab = jbeg, jend
661 <                j = list(jnab)
764 <                              
765 <                call get_interatomic_vector(q_group_Row(:,i), &
766 <                     q_group_Col(:,j), d_grp, rgrpsq)
767 <                
768 <                vab = 0.0d0
769 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
770 <                     in_switching_region)
771 <                
772 <                do ia = groupStart(i), groupStart(i+1)-1
773 <                   atom1 = groupList(ia)
774 <                  
775 <                   do jb = groupStart(j), groupStart(j+1)-1
776 <                      atom2 = groupList(jb)                        
777 <                      
778 <                      call get_interatomic_vector(q_Row(:,atom1), &
779 <                           q_Col(:,atom2), d_atm, ratmsq)
780 <                                            
781 <                      call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
782 <                           do_pot, do_stress, &
783 <                           u_l, A, f, t, pot_local, vpair)
784 <                      
785 <                      vab = vab + vpair
786 <                                            
787 <                   enddo
788 <                enddo
789 <                
790 <                if (in_switching_region) then
791 <                   swderiv = vab*dswdr/rgrp
792 <  
793 <                   do ia=groupStart(i), groupStart(i+1)-1
794 <                      atom1=groupList(ia)
795 <                      mf = mfact(atom1)                  
796 <                      f_Row(1,atom1) = f_Row(1,atom1) - swderiv*d_grp(1)*mf
797 <                      f_Row(2,atom1) = f_Row(2,atom1) - swderiv*d_grp(2)*mf
798 <                      f_Row(3,atom1) = f_Row(3,atom1) - swderiv*d_grp(3)*mf
799 <                   enddo
800 <                  
801 <                   do jb=groupStart(j), groupStart(j+1)-1
802 <                      atom2=groupList(jb)
803 <                      mf = mfact(atom2)
804 <                      f_Col(1,atom2) = f_Col(1,atom2) + swderiv*d_grp(1)*mf
805 <                      f_Col(2,atom2) = f_Col(2,atom2) + swderiv*d_grp(2)*mf
806 <                      f_Col(3,atom2) = f_Col(3,atom2) + swderiv*d_grp(3)*mf
807 <                   enddo
808 <                endif
809 <
810 <             enddo
651 >       if (update_nlist) then
652 > #ifdef IS_MPI
653 >          point(nrow_group + 1) = nlist + 1
654 > #else
655 >          point(nGroup) = nlist + 1
656 > #endif
657 >          if (loop .eq. PREPAIR_LOOP) then
658 >             ! we just did the neighbor list update on the first
659 >             ! pass, so we don't need to do it
660 >             ! again on the second pass
661 >             update_nlist = .false.                              
662            endif
663 <       enddo
813 <    endif
814 <    
815 < #else
816 <
817 <    if (update_nlist) then
818 <      
819 <       !! save current configuration, construct neighbor list,
820 <       !! and calculate forces
821 <      
822 <       call saveNeighborList(nGroup, q_group)
823 <      
824 <       neighborListSize = size(list)
825 <       nlist = 0      
826 <      
827 <       do i = 1, nGroup-1
828 <          point(i) = nlist + 1
829 <          
830 <          do j = i+1, nGroup
663 >       endif
664              
665 <             call get_interatomic_vector(q_group(:,i), &
666 <                  q_group(:,j), d_grp, rgrpsq)
667 <            
835 <             if (rgrpsq < rlistsq) then
836 <                nlist = nlist + 1
837 <                
838 <                if (nlist > neighborListSize) then
839 <                   call expandNeighborList(nGroup, listerror)
840 <                   if (listerror /= 0) then
841 <                      error = -1
842 <                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
843 <                      return
844 <                   end if
845 <                   neighborListSize = size(list)
846 <                endif
847 <                
848 <                list(nlist) = j
849 <
850 <                vab = 0.0d0
851 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
852 <                     in_switching_region)
853 <                
854 <                do ia = groupStart(i), groupStart(i+1)-1
855 <                   atom1 = groupList(ia)
856 <                  
857 <                   inner: do jb = groupStart(j), groupStart(j+1)-1
858 <                      atom2 = groupList(jb)
859 <                      
860 <                      if (skipThisPair(atom1, atom2)) cycle inner
861 <                      
862 <                      call get_interatomic_vector(q(:,atom1), &
863 <                           q(:,atom2), d_atm, ratmsq)
864 <                      
865 <                      call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
866 <                           do_pot, do_stress, &
867 <                           u_l, A, f, t, pot, vpair)
868 <                      
869 <                      vab = vab + vpair
870 <                      
871 <                   enddo inner
872 <                enddo
873 <
874 <                if (in_switching_region) then
875 <                   swderiv = vab*dswdr/rgrp
876 <                   do ia=groupStart(i), groupStart(i+1)-1
877 <                      atom1=groupList(ia)
878 <                      mf = mfact(atom1)                  
879 <                      f(1,atom1) = f(1,atom1) - swderiv*d_grp(1)*mf
880 <                      f(2,atom1) = f(2,atom1) - swderiv*d_grp(2)*mf
881 <                      f(3,atom1) = f(3,atom1) - swderiv*d_grp(3)*mf
882 <                   enddo
883 <                  
884 <                   do jb=groupStart(j), groupStart(j+1)-1
885 <                      atom2=groupList(jb)
886 <                      mf = mfact(atom2)
887 <                      f(1,atom2) = f(1,atom2) + swderiv*d_grp(1)*mf
888 <                      f(2,atom2) = f(2,atom2) + swderiv*d_grp(2)*mf
889 <                      f(3,atom2) = f(3,atom2) + swderiv*d_grp(3)*mf
890 <                   enddo
891 <                endif
892 <
893 <             end if
894 <          enddo
895 <       enddo
896 <       point(nGroup) = nlist + 1          
665 >       if (loop .eq. PREPAIR_LOOP) then
666 >          call do_preforce(nlocal, pot)
667 >       endif
668        
669 <    else  !! (of update_check)
899 <      
900 <       ! use the list to find the neighbors
901 <       do i = 1, nGroup-1
902 <          JBEG = POINT(i)
903 <          JEND = POINT(i+1) - 1
904 <          ! check that group i has neighbors
905 <          if (jbeg .le. jend) then
906 <            
907 <             do jnab = jbeg, jend
908 <                j = list(jnab)
909 <
910 <                call get_interatomic_vector(q_group(:,i), &
911 <                     q_group(:,j), d_grp, rgrpsq)
912 <
913 <                vab = 0.0d0                    
914 <                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
915 <                     in_switching_region)
916 <                                
917 <                do ia = groupStart(i), groupStart(i+1)-1
918 <                   atom1 = groupList(ia)
919 <                  
920 <                   do jb = groupStart(j), groupStart(j+1)-1
921 <                      atom2 = groupList(jb)                        
922 <                      
923 <                      call get_interatomic_vector(q(:,atom1), &
924 <                           q(:,atom2), d_atm, ratmsq)
925 <                      
926 <                      call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
927 <                           do_pot, do_stress, &
928 <                           u_l, A, f, t, pot, vpair)
929 <                      
930 <                      vab = vab + vpair
931 <                      
932 <                   enddo
933 <                enddo
934 <                
935 <                if (in_switching_region) then
936 <                   swderiv = vab*dswdr/rgrp
937 <                  
938 <                   do ia=groupStart(i), groupStart(i+1)-1
939 <                      atom1=groupList(ia)
940 <                      mf = mfact(atom1)                  
941 <                      f(1,atom1) = f(1,atom1) - swderiv*d_grp(1)*mf
942 <                      f(2,atom1) = f(2,atom1) - swderiv*d_grp(2)*mf
943 <                      f(3,atom1) = f(3,atom1) - swderiv*d_grp(3)*mf
944 <                   enddo
945 <                  
946 <                   do jb=groupStart(j), groupStart(j+1)-1
947 <                      atom2=groupList(jb)
948 <                      mf = mfact(atom2)
949 <                      f(1,atom2) = f(1,atom2) + swderiv*d_grp(1)*mf
950 <                      f(2,atom2) = f(2,atom2) + swderiv*d_grp(2)*mf
951 <                      f(3,atom2) = f(3,atom2) + swderiv*d_grp(3)*mf
952 <                   enddo
953 <                endif
954 <             enddo
955 <          endif
956 <       enddo
957 <    endif
958 <        
959 < #endif
669 >    enddo
670      
961 –    ! phew, done with main loop.
962 –    
671      !! Do timing
672   #ifdef PROFILE
673      call cpu_time(forceTimeFinal)
# Line 1081 | Line 789 | contains
789      endif
790      
791   #endif
792 <    
1085 <    
792 >      
793    end subroutine do_force_loop
1087 –
794    
795 <  subroutine do_pair(i, j, rijsq, d, sw, do_pot, do_stress, &
796 <       u_l, A, f, t, pot, vpair)
795 >  subroutine do_pair(i, j, rijsq, d, sw, do_pot, &
796 >       u_l, A, f, t, pot, vpair, fpair)
797  
798      real( kind = dp ) :: pot, vpair, sw
799 +    real( kind = dp ), dimension(3) :: fpair
800      real( kind = dp ), dimension(nLocal)   :: mfact
801      real( kind = dp ), dimension(3,nLocal) :: u_l
802      real( kind = dp ), dimension(9,nLocal) :: A
803      real( kind = dp ), dimension(3,nLocal) :: f
804      real( kind = dp ), dimension(3,nLocal) :: t
805  
806 <    logical, intent(inout) :: do_pot, do_stress
806 >    logical, intent(inout) :: do_pot
807      integer, intent(in) :: i, j
808      real ( kind = dp ), intent(inout) :: rijsq
809      real ( kind = dp )                :: r
# Line 1104 | Line 811 | contains
811      integer :: me_i, me_j
812  
813      r = sqrt(rijsq)
814 +    vpair = 0.0d0
815 +    fpair(1:3) = 0.0d0
816  
817   #ifdef IS_MPI
818      if (tagRow(i) .eq. tagColumn(j)) then
# Line 1119 | Line 828 | contains
828      if (FF_uses_LJ .and. SIM_uses_LJ) then
829        
830         if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) then
831 <          call do_lj_pair(i, j, d, r, rijsq, sw, vpair, pot, f, do_pot, &
832 <               do_stress)
831 >          !write(*,*) 'calling lj with'
832 >          !write(*,*) i, j, r, rijsq
833 >          !write(*,'(3es12.3)') d(1), d(2), d(3)
834 >          !write(*,'(3es12.3)') sw, vpair, pot
835 >          !write(*,*)
836 >
837 >          call do_lj_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
838         endif
839        
840      endif
# Line 1128 | Line 842 | contains
842      if (FF_uses_charges .and. SIM_uses_charges) then
843        
844         if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then
845 <          call do_charge_pair(i, j, d, r, rijsq, sw, vpair, pot, f, do_pot, &
1132 <               do_stress)
845 >          call do_charge_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, do_pot)
846         endif
847        
848      endif
# Line 1137 | Line 850 | contains
850      if (FF_uses_dipoles .and. SIM_uses_dipoles) then
851        
852         if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
853 <          call do_dipole_pair(i, j, d, r, rijsq, sw, vpair, pot, u_l, f, t, &
854 <               do_pot, do_stress)
853 >          call do_dipole_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, u_l, f, t, &
854 >               do_pot)
855            if (FF_uses_RF .and. SIM_uses_RF) then
856               call accumulate_rf(i, j, r, u_l, sw)
857 <             call rf_correct_forces(i, j, d, r, u_l, sw, f, do_stress)
857 >             call rf_correct_forces(i, j, d, r, u_l, sw, f, fpair)
858            endif          
859         endif
860  
# Line 1150 | Line 863 | contains
863      if (FF_uses_Sticky .and. SIM_uses_sticky) then
864  
865         if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
866 <          call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, pot, A, f, t, &
867 <               do_pot, do_stress)
866 >          call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, A, f, t, &
867 >               do_pot)
868         endif
869  
870      endif
# Line 1160 | Line 873 | contains
873      if (FF_uses_GB .and. SIM_uses_GB) then
874        
875         if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
876 <          call do_gb_pair(i, j, d, r, rijsq, sw, vpair, pot, u_l, f, t, &
877 <               do_pot, do_stress)          
876 >          call do_gb_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, u_l, f, t, &
877 >               do_pot)
878         endif
879  
880      endif
# Line 1169 | Line 882 | contains
882      if (FF_uses_EAM .and. SIM_uses_EAM) then
883        
884         if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
885 <          call do_eam_pair(i, j, d, r, rijsq, sw, vpair, pot, f, &
886 <               do_pot, do_stress)
885 >          call do_eam_pair(i, j, d, r, rijsq, sw, vpair, fpair, pot, f, &
886 >               do_pot)
887         endif
888        
889      endif
890      
891    end subroutine do_pair
892  
893 <  subroutine do_prepair(i, j, rijsq, d, rcijsq, dc, &
893 >  subroutine do_prepair(i, j, rijsq, d, sw, rcijsq, dc, &
894         do_pot, do_stress, u_l, A, f, t, pot)
895 <   real( kind = dp ) :: pot
895 >
896 >   real( kind = dp ) :: pot, sw
897     real( kind = dp ), dimension(3,nLocal) :: u_l
898     real (kind=dp), dimension(9,nLocal) :: A
899     real (kind=dp), dimension(3,nLocal) :: f
# Line 1206 | Line 920 | contains
920  
921   #ifdef IS_MPI
922     if (tagRow(i) .eq. tagColumn(j)) then
923 <      write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
923 >      write(0,*) 'do_prepair is doing', i , j, tagRow(i), tagColumn(j)
924     endif
925    
926     me_i = atid_row(i)
# Line 1218 | Line 932 | contains
932     me_j = atid(j)
933    
934   #endif
935 <    
935 >  
936     if (FF_uses_EAM .and. SIM_uses_EAM) then
937 <
937 >      
938        if (PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) &
939             call calc_EAM_prepair_rho(i, j, d, r, rijsq )
940 <
940 >      
941     endif
942    
943   end subroutine do_prepair
944 <
945 <
1232 <
1233 <
944 >
945 >
946   subroutine do_preforce(nlocal,pot)
947     integer :: nlocal
948     real( kind = dp ) :: pot
# Line 1409 | Line 1121 | contains
1121        endif
1122     enddo
1123    
1124 <   do i = 1, nExcludes_local
1125 <      if (excludesLocal(1,i) == unique_id_1) then
1126 <         if (excludesLocal(2,i) == unique_id_2) then
1127 <            skip_it = .true.
1416 <            return
1417 <         endif
1418 <      else
1419 <         if (excludesLocal(1,i) == unique_id_2) then
1420 <            if (excludesLocal(2,i) == unique_id_1) then
1421 <               skip_it = .true.
1422 <               return
1423 <            endif
1424 <         endif
1124 >   do i = 1, nSkipsForAtom(unique_id_1)
1125 >      if (skipsForAtom(unique_id_1, i) .eq. unique_id_2) then
1126 >         skip_it = .true.
1127 >         return
1128        endif
1129     end do
1130    
# Line 1452 | Line 1155 | contains
1155   #endif
1156  
1157   !! This cleans componets of force arrays belonging only to fortran
1158 +
1159 + subroutine add_stress_tensor(dpair, fpair)
1160 +  
1161 +   real( kind = dp ), dimension(3), intent(in) :: dpair, fpair
1162 +  
1163 +   ! because the d vector is the rj - ri vector, and
1164 +   ! because fx, fy, fz are the force on atom i, we need a
1165 +   ! negative sign here:  
1166 +  
1167 +   tau_Temp(1) = tau_Temp(1) - dpair(1) * fpair(1)
1168 +   tau_Temp(2) = tau_Temp(2) - dpair(1) * fpair(2)
1169 +   tau_Temp(3) = tau_Temp(3) - dpair(1) * fpair(3)
1170 +   tau_Temp(4) = tau_Temp(4) - dpair(2) * fpair(1)
1171 +   tau_Temp(5) = tau_Temp(5) - dpair(2) * fpair(2)
1172 +   tau_Temp(6) = tau_Temp(6) - dpair(2) * fpair(3)
1173 +   tau_Temp(7) = tau_Temp(7) - dpair(3) * fpair(1)
1174 +   tau_Temp(8) = tau_Temp(8) - dpair(3) * fpair(2)
1175 +   tau_Temp(9) = tau_Temp(9) - dpair(3) * fpair(3)
1176 +  
1177 +   !write(*,'(6es12.3)')  fpair(1:3), tau_Temp(1), tau_Temp(5), tau_temp(9)
1178 +   virial_Temp = virial_Temp + &
1179 +        (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))
1180 +  
1181 + end subroutine add_stress_tensor
1182  
1183   end module do_Forces
1184 +

Diff Legend

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