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 1144 by tim, Sat May 1 18:52:38 2004 UTC vs.
Revision 1150 by gezelter, Fri May 7 21:35:05 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.53 2004-05-01 18:52:38 tim Exp $, $Date: 2004-05-01 18:52:38 $, $Name: not supported by cvs2svn $, $Revision: 1.53 $
7 > !! @version $Id: do_Forces.F90,v 1.54 2004-05-07 21:35:04 gezelter Exp $, $Date: 2004-05-07 21:35:04 $, $Name: not supported by cvs2svn $, $Revision: 1.54 $
8  
9   module do_Forces
10    use force_globals
11    use simulation
12    use definitions
13    use atype_module
14 +  use switcheroo
15    use neighborLists  
16    use lj
17    use sticky_pair
# Line 30 | Line 31 | module do_Forces
31  
32   #define __FORTRAN90
33   #include "fForceField.h"
34 + #include "fSwitchingFunction.h"
35  
36    logical, save :: haveRlist = .false.
37    logical, save :: haveNeighborList = .false.
# Line 62 | Line 64 | module do_Forces
64    public :: init_FF
65    public :: do_force_loop
66    public :: setRlistDF
65
67  
68   #ifdef PROFILE
69    public :: getforcetime
# Line 364 | Line 365 | contains
365         endif
366         haveNeighborList = .true.
367      endif
368 +
369      
370 +    
371    end subroutine init_FF
372    
373  
374    !! Does force loop over i,j pairs. Calls do_pair to calculates forces.
375    !------------------------------------------------------------->
376 <  subroutine do_force_loop(q, qcom, A, u_l, f, t, tau, pot, &
376 >  subroutine do_force_loop(q, q_group, A, u_l, f, t, tau, pot, &
377         do_pot_c, do_stress_c, error)
378      !! Position array provided by C, dimensioned by getNlocal
379      real ( kind = dp ), dimension(3,nLocal) :: q
380      !! molecular center-of-mass position array
381 <    real ( kind = dp ), dimension(3,nLocal) :: qcom
381 >    real ( kind = dp ), dimension(3,nGroup) :: q_group
382      !! Rotation Matrix for each long range particle in simulation.
383      real( kind = dp), dimension(9,nLocal) :: A    
384      !! Unit vectors for dipoles (lab frame)
# Line 391 | Line 394 | contains
394      logical ( kind = 2) :: do_pot_c, do_stress_c
395      logical :: do_pot
396      logical :: do_stress
397 +    logical :: in_switching_region
398   #ifdef IS_MPI
399      real( kind = DP ) :: pot_local
400      integer :: nrow
401      integer :: ncol
402      integer :: nprocs
403 +    integer :: nrow_group
404 +    integer :: ncol_group
405   #endif
406      integer :: natoms    
407      logical :: update_nlist  
408      integer :: i, j, jbeg, jend, jnab
409 +    integer :: ia, jb, atom1, atom2
410      integer :: nlist
411 <    real( kind = DP ) ::  rijsq, rcijsq
412 <    real(kind=dp),dimension(3) :: d, dc
411 >    real( kind = DP ) :: ratmsq, rgrpsq, rgrp, vpair, vab
412 >    real( kind = DP ) :: sw, dswdr, swderiv, mf
413 >    real(kind=dp),dimension(3) :: d_atm, d_grp
414      real(kind=dp) :: rfpot, mu_i, virial
415      integer :: me_i, me_j
416      logical :: is_dp_i
# Line 419 | Line 427 | contains
427      pot_local = 0.0_dp
428      nrow   = getNrow(plan_row)
429      ncol   = getNcol(plan_col)
430 +    nrow_group   = getNrowGroup(plan_row)
431 +    ncol_group   = getNcolGroup(plan_col)
432   #else
433      natoms = nlocal
434   #endif
# Line 438 | Line 448 | contains
448      
449   #ifdef IS_MPI    
450      
451 <    call gather(q,q_Row,plan_row3d)
452 <    call gather(q,q_Col,plan_col3d)
453 <    
454 <    if (SIM_uses_molecular_cutoffs) then
455 <       call gather(qcom,qcom_Row,plan_row3d)
456 <       call gather(qcom,qcom_Col,plan_col3d)
447 <    endif
448 <    
451 >    call gather(q, q_Row, plan_row3d)
452 >    call gather(q, q_Col, plan_col3d)
453 >
454 >    call gather(q_group, q_group_Row, plan_row_Group_3d)
455 >    call gather(q_group, q_group_Col, plan_col_Group_3d)
456 >        
457      if (FF_UsesDirectionalAtoms() .and. SIM_uses_directional_atoms) then
458         call gather(u_l,u_l_Row,plan_row3d)
459         call gather(u_l,u_l_Col,plan_col3d)
# Line 464 | Line 472 | contains
472      
473      if (FF_RequiresPrepairCalc() .and. SIM_requires_prepair_calc) then
474         !! See if we need to update neighbor lists
475 <       if (SIM_uses_molecular_cutoffs) then
476 <          call checkNeighborList(nlocal, qcom, listSkin, update_nlist)
477 <       else
470 <          call checkNeighborList(nlocal, q, listSkin, update_nlist)  
471 <       endif
475 >
476 >       call checkNeighborList(nGroup, q_group, listSkin, update_nlist)
477 >
478         !! if_mpi_gather_stuff_for_prepair
479         !! do_prepair_loop_if_needed
480         !! if_mpi_scatter_stuff_from_prepair
# Line 481 | Line 487 | contains
487            
488            !! save current configuration, construct neighbor list,
489            !! and calculate forces
490 <          if (SIM_uses_molecular_cutoffs) then
491 <             call saveNeighborList(nlocal, qcom)
486 <          else
487 <             call saveNeighborList(nlocal, q)
488 <          endif
490 >
491 >          call saveNeighborList(nGroup, q_group)
492            
493            neighborListSize = size(list)
494            nlist = 0      
495            
496 <          do i = 1, nrow
496 >          do i = 1, nrow_group
497               point(i) = nlist + 1
498              
499 <             prepair_inner: do j = 1, ncol
499 >             do j = 1, ncol_group
500                  
501 <                if (skipThisPair(i,j)) cycle prepair_inner
501 >                call get_interatomic_vector(q_group_Row(:,i), &
502 >                     q_group_Col(:,j), d_grp, rgrpsq)
503                  
504 <                if (SIM_uses_molecular_cutoffs) then
501 <                   call get_interatomic_vector(qcom_Row(:,i), qcom_Col(:,j), &
502 <                        dc, rcijsq)
503 <                else
504 <                   call get_interatomic_vector(q_Row(:,i), q_Col(:,j), &
505 <                        dc, rcijsq)
506 <                endif
507 <                
508 <                if (rcijsq < rlistsq) then            
509 <                  
504 >                if (rgrpsq < rlistsq) then
505                     nlist = nlist + 1
506                    
507                     if (nlist > neighborListSize) then
508 <                      call expandNeighborList(nlocal, listerror)
508 >                      call expandNeighborList(nGroup, listerror)
509                        if (listerror /= 0) then
510                           error = -1
511                           write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
# Line 521 | Line 516 | contains
516                    
517                     list(nlist) = j
518                    
519 <                   if (SIM_uses_molecular_cutoffs) then
520 <                      ! since the simulation distances were in molecular COMs:
521 <                      call get_interatomic_vector(q_Row(:,i), q_Col(:,j), &
522 <                           d, rijsq)
523 <                   else
524 <                      d(1:3) = dc(1:3)
525 <                      rijsq = rcijsq
526 <                   endif
527 <                  
528 <                   call do_prepair(i, j, rijsq, d, rcijsq, dc, &
529 <                        do_pot, do_stress, u_l, A, f, t, pot_local)
530 <                endif
531 <             enddo prepair_inner
519 >                   do ia = groupStart(i), groupStart(i+1)-1
520 >                      atom1 = groupList(ia)
521 >
522 >                      prepair_inner: do jb = groupStart(j), groupStart(j+1)-1
523 >                         atom2 = groupList(jb)
524 >                        
525 >                         if (skipThisPair(atom1, atom2)) cycle prepair_inner
526 >                        
527 >                         call get_interatomic_vector(q_Row(:,atom1), &
528 >                              q_Col(:,atom2), d_atm, ratmsq)
529 >                        
530 >                         call do_prepair(atom1, atom2, ratmsq, d_atm, &
531 >                              rgrpsq, d_grp, do_pot, do_stress, &
532 >                              u_l, A, f, t, pot_local)
533 >
534 >                      enddo prepair_inner
535 >                   enddo
536 >                end if
537 >             enddo
538            enddo
539 <          
540 <          point(nrow + 1) = nlist + 1
540 <          
539 >          point(nrow_group + 1) = nlist + 1          
540 >                          
541         else  !! (of update_check)
542            
543            ! use the list to find the neighbors
544 <          do i = 1, nrow
544 >          do i = 1, nrow_group
545               JBEG = POINT(i)
546               JEND = POINT(i+1) - 1
547 <             ! check thiat molecule i has neighbors
547 >             ! check that group i has neighbors
548               if (jbeg .le. jend) then
549                  
550                  do jnab = jbeg, jend
551                     j = list(jnab)
552                                      
553                   call get_interatomic_vector(q_Row(:,i), q_Col(:,j), &
554                        d, rijsq)
555                   if (SIM_uses_molecular_cutoffs) then
556                      call get_interatomic_vector(qcom_Row(:,i),qcom_Col(:,j),&
557                           dc, rcijsq)
558                   else
559                      dc(1:3) = d(1:3)
560                      rcijsq = rijsq
561                   endif
552                    
553 <                   call do_prepair(i, j, rijsq, d, rcijsq, dc, &
554 <                        do_pot, do_stress, u_l, A, f, t, pot_local)
555 <                  
553 >                   do ia = groupStart(i), groupStart(i+1)-1
554 >                      atom1 = groupList(ia)
555 >
556 >                      do jb = groupStart(j), groupStart(j+1)-1
557 >                         atom2 = groupList(jb)                        
558 >                        
559 >                         call get_interatomic_vector(q_Row(:,atom1), &
560 >                              q_Col(:,atom2), d_atm, ratmsq)
561 >                        
562 >                         call do_prepair(atom1, atom2, ratmsq, d_atm, &
563 >                              rgrpsq, d_grp, do_pot, do_stress, &
564 >                              u_l, A, f, t, pot_local)
565 >                        
566 >                      enddo
567 >                   enddo
568                  enddo
569               endif
570            enddo
571         endif
572        
573   #else
574 <      
574 >
575         if (update_nlist) then
576            
577 <          ! save current configuration, contruct neighbor list,
578 <          ! and calculate forces
577 <          call saveNeighborList(natoms, q)
577 >          !! save current configuration, construct neighbor list,
578 >          !! and calculate forces
579            
580 +          call saveNeighborList(nGroup, q_group)
581 +          
582            neighborListSize = size(list)
583 +          nlist = 0      
584            
585 <          nlist = 0
582 <          
583 <          do i = 1, natoms-1
585 >          do i = 1, nGroup-1
586               point(i) = nlist + 1
587              
588 <             prepair_inner: do j = i+1, natoms
588 >             do j = i+1, nGroup
589                  
590 <                if (skipThisPair(i,j))  cycle prepair_inner
590 >                call get_interatomic_vector(q_group(:,i), &
591 >                     q_group(:,j), d_grp, rgrpsq)
592                  
593 <                if (SIM_uses_molecular_cutoffs) then
591 <                   call get_interatomic_vector(qcom(:,i), qcom(:,j), &
592 <                        dc, rcijsq)
593 <                else
594 <                   call get_interatomic_vector(q(:,i), q(:,j), dc, rcijsq)
595 <                endif
596 <                
597 <                if (rcijsq < rlistsq) then
598 <                  
593 >                if (rgrpsq < rlistsq) then
594                     nlist = nlist + 1
595                    
596                     if (nlist > neighborListSize) then
597 <                      call expandNeighborList(natoms, listerror)
597 >                      call expandNeighborList(nGroup, listerror)
598                        if (listerror /= 0) then
599                           error = -1
600                           write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
# Line 610 | Line 605 | contains
605                    
606                     list(nlist) = j
607                    
608 <                  
609 <                   if (SIM_uses_molecular_cutoffs) then
610 <                      ! since the simulation distances were in molecular COMs:
611 <                      call get_interatomic_vector(q(:,i), q(:,j), &
612 <                           d, rijsq)
613 <                   else
614 <                      dc(1:3) = d(1:3)
615 <                      rcijsq = rijsq
616 <                   endif
617 <                  
618 <                   call do_prepair(i, j, rijsq, d, rcijsq, dc, &
619 <                        do_pot, do_stress, u_l, A, f, t, pot)
620 <                  
621 <                endif
622 <             enddo prepair_inner
608 >                   do ia = groupStart(i), groupStart(i+1)-1
609 >                      atom1 = groupList(ia)
610 >                      
611 >                      prepair_inner: do jb = groupStart(j), groupStart(j+1)-1
612 >                         atom2 = groupList(jb)
613 >                        
614 >                         if (skipThisPair(atom1, atom2)) cycle prepair_inner
615 >                        
616 >                         call get_interatomic_vector(q(:,atom1), &
617 >                              q(:,atom2), d_atm, ratmsq)
618 >                        
619 >                         call do_prepair(atom1, atom2, ratmsq, d_atm, &
620 >                              rgrpsq, d_grp, do_pot, do_stress, &
621 >                              u_l, A, f, t, pot)
622 >                        
623 >                      enddo prepair_inner
624 >                   enddo
625 >                end if
626 >             enddo
627            enddo
628 +          point(nGroup) = nlist + 1          
629 +                          
630 +       else  !! (of update_check)
631            
630          point(natoms) = nlist + 1
631          
632       else !! (update)
633          
632            ! use the list to find the neighbors
633 <          do i = 1, natoms-1
633 >          do i = 1, nGroup-1
634               JBEG = POINT(i)
635               JEND = POINT(i+1) - 1
636 <             ! check thiat molecule i has neighbors
636 >             ! check that group i has neighbors
637               if (jbeg .le. jend) then
638                  
639                  do jnab = jbeg, jend
640                     j = list(jnab)
641                    
642 <                   call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
643 <                   if (SIM_uses_molecular_cutoffs) then
644 <                      call get_interatomic_vector(qcom(:,i), qcom(:,j), &
645 <                           dc, rcijsq)
646 <                   else
647 <                      dc(1:3) = d(1:3)
648 <                      rcijsq = rijsq
649 <                   endif
650 <                  
651 <                   call do_prepair(i, j, rijsq, d, rcijsq, dc, &
652 <                        do_pot, do_stress, u_l, A, f, t, pot)
653 <                  
642 >                   do ia = groupStart(i), groupStart(i+1)-1
643 >                      atom1 = groupList(ia)
644 >
645 >                      do jb = groupStart(j), groupStart(j+1)-1
646 >                         atom2 = groupList(jb)                        
647 >                        
648 >                         call get_interatomic_vector(q(:,atom1), &
649 >                              q(:,atom2), d_atm, ratmsq)
650 >                        
651 >                         call do_prepair(atom1, atom2, ratmsq, d_atm, &
652 >                              rgrpsq, d_grp, do_pot, do_stress, &
653 >                              u_l, A, f, t, pot)
654 >                        
655 >                      enddo
656 >                   enddo
657                  enddo
658               endif
659            enddo
660         endif
661 +      
662   #endif
663 +      
664         !! Do rest of preforce calculations
665         !! do necessary preforce calculations  
666         call do_preforce(nlocal,pot)
# Line 665 | Line 668 | contains
668         update_nlist = .false.
669      else
670         !! See if we need to update neighbor lists for non pre-pair
671 <       call checkNeighborList(nlocal, q, listSkin, update_nlist)  
671 >       call checkNeighborList(nGroup, q_group, listSkin, update_nlist)  
672      endif
673      
674      !---------------------------------MAIN Pair LOOP->>>>>>>>>>>>>
# Line 673 | Line 676 | contains
676   #ifdef IS_MPI
677      
678      if (update_nlist) then
679 +      
680         !! save current configuration, construct neighbor list,
681         !! and calculate forces
678       call saveNeighborList(nlocal, q)
682        
683 +       call saveNeighborList(nGroup, q_group)
684 +      
685         neighborListSize = size(list)
686         nlist = 0      
687        
688 <       do i = 1, nrow
684 <          
688 >       do i = 1, nrow_group
689            point(i) = nlist + 1
690            
691 <          inner: do j = 1, ncol
691 >          do j = 1, ncol_group
692              
693 <             if (skipThisPair(i,j)) cycle inner
693 >             call get_interatomic_vector(q_group_Row(:,i), &
694 >                  q_group_Col(:,j), d_grp, rgrpsq)
695              
696 <             if (SIM_uses_molecular_cutoffs) then
692 <                call get_interatomic_vector(qcom_Row(:,i), qcom_Col(:,j), &
693 <                     dc, rcijsq)
694 <             else
695 <                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), &
696 <                     dc, rcijsq)
697 <             endif
698 <            
699 <             if (rcijsq < rlistsq) then            
700 <                
696 >             if (rgrpsq < rlistsq) then
697                  nlist = nlist + 1
698                  
699                  if (nlist > neighborListSize) then
700 <                   call expandNeighborList(nlocal, listerror)
700 >                   call expandNeighborList(nGroup, listerror)
701                     if (listerror /= 0) then
702                        error = -1
703                        write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
# Line 712 | Line 708 | contains
708                  
709                  list(nlist) = j
710                  
711 <                if (SIM_uses_molecular_cutoffs) then
712 <                   call get_interatomic_vector(q_Row(:,i), q_Col(:,j), &
713 <                        d, rijsq)
718 <                else
719 <                   d(1:3) = dc(1:3)
720 <                   rijsq = rcijsq
721 <                endif
711 >                vab = 0.0d0
712 >                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
713 >                     in_switching_region)
714                  
715 <                call do_pair(i, j, rijsq, d, rcijsq, dc, mfact, &
716 <                     do_pot, do_stress, u_l, A, f, t, pot_local)
717 <                
718 <             endif
719 <          enddo inner
715 >                do ia = groupStart(i), groupStart(i+1)-1
716 >                   atom1 = groupList(ia)
717 >                  
718 >                   inner: do jb = groupStart(j), groupStart(j+1)-1
719 >                      atom2 = groupList(jb)
720 >                      
721 >                      if (skipThisPair(atom1, atom2)) cycle inner
722 >                      
723 >                      call get_interatomic_vector(q_Row(:,atom1), &
724 >                           q_Col(:,atom2), d_atm, ratmsq)
725 >                      
726 >                      call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
727 >                           do_pot, do_stress, &
728 >                           u_l, A, f, t, pot_local, vpair)
729 >                      
730 >                      vab = vab + vpair
731 >                   enddo inner
732 >                enddo
733 >
734 >                if (in_switching_region) then
735 >                   swderiv = vab*dswdr/rgrp
736 >                  
737 >                   do ia=groupStart(i), groupStart(i+1)-1
738 >                      atom1=groupList(ia)
739 >                      mf = mfact(atom1)                  
740 >                      f_Row(1,atom1) = f_Row(1,atom1) - swderiv*d_grp(1)*mf
741 >                      f_Row(2,atom1) = f_Row(2,atom1) - swderiv*d_grp(2)*mf
742 >                      f_Row(3,atom1) = f_Row(3,atom1) - swderiv*d_grp(3)*mf
743 >                   enddo
744 >                  
745 >                   do jb=groupStart(j), groupStart(j+1)-1
746 >                      atom2=groupList(jb)
747 >                      mf = mfact(atom2)
748 >                      f_Col(1,atom2) = f_Col(1,atom2) + swderiv*d_grp(1)*mf
749 >                      f_Col(2,atom2) = f_Col(2,atom2) + swderiv*d_grp(2)*mf
750 >                      f_Col(3,atom2) = f_Col(3,atom2) + swderiv*d_grp(3)*mf
751 >                   enddo                        
752 >                endif
753 >
754 >             end if
755 >          enddo
756         enddo
757 +       point(nrow_group + 1) = nlist + 1          
758        
730       point(nrow + 1) = nlist + 1
731      
759      else  !! (of update_check)
760        
761         ! use the list to find the neighbors
762 <       do i = 1, nrow
762 >       do i = 1, nrow_group
763            JBEG = POINT(i)
764            JEND = POINT(i+1) - 1
765 <          ! check thiat molecule i has neighbors
765 >          ! check that group i has neighbors
766            if (jbeg .le. jend) then
767              
768               do jnab = jbeg, jend
769                  j = list(jnab)
770 +                              
771 +                call get_interatomic_vector(q_group_Row(:,i), &
772 +                     q_group_Col(:,j), d_grp, rgrpsq)
773                  
774 <                call get_interatomic_vector(q_Row(:,i), q_Col(:,j), d, rijsq)
775 <                if (SIM_uses_molecular_cutoffs) then
776 <                   call get_interatomic_vector(qcom_Row(:,i), qcom_Col(:,j), &
747 <                        dc, rcijsq)
748 <                else
749 <                   dc(1:3) = d(1:3)
750 <                   rcijsq = rijsq
751 <                endif
774 >                vab = 0.0d0
775 >                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
776 >                     in_switching_region)
777                  
778 <                call do_pair(i, j, rijsq, d, rcijsq, dc, mfact, &
779 <                     do_pot, do_stress, u_l, A, f, t, pot_local)
778 >                do ia = groupStart(i), groupStart(i+1)-1
779 >                   atom1 = groupList(ia)
780 >                  
781 >                   do jb = groupStart(j), groupStart(j+1)-1
782 >                      atom2 = groupList(jb)                        
783 >                      
784 >                      call get_interatomic_vector(q_Row(:,atom1), &
785 >                           q_Col(:,atom2), d_atm, ratmsq)
786 >                                            
787 >                      call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
788 >                           do_pot, do_stress, &
789 >                           u_l, A, f, t, pot_local, vpair)
790 >                      
791 >                      vab = vab + vpair
792 >                                            
793 >                   enddo
794 >                enddo
795                  
796 +                if (in_switching_region) then
797 +                   swderiv = vab*dswdr/rgrp
798 +  
799 +                   do ia=groupStart(i), groupStart(i+1)-1
800 +                      atom1=groupList(ia)
801 +                      mf = mfact(atom1)                  
802 +                      f_Row(1,atom1) = f_Row(1,atom1) - swderiv*d_grp(1)*mf
803 +                      f_Row(2,atom1) = f_Row(2,atom1) - swderiv*d_grp(2)*mf
804 +                      f_Row(3,atom1) = f_Row(3,atom1) - swderiv*d_grp(3)*mf
805 +                   enddo
806 +                  
807 +                   do jb=groupStart(j), groupStart(j+1)-1
808 +                      atom2=groupList(jb)
809 +                      mf = mfact(atom2)
810 +                      f_Col(1,atom2) = f_Col(1,atom2) + swderiv*d_grp(1)*mf
811 +                      f_Col(2,atom2) = f_Col(2,atom2) + swderiv*d_grp(2)*mf
812 +                      f_Col(3,atom2) = f_Col(3,atom2) + swderiv*d_grp(3)*mf
813 +                   enddo
814 +                endif
815 +
816               enddo
817            endif
818         enddo
819      endif
820      
821   #else
822 <    
822 >
823      if (update_nlist) then
824        
825 <       ! save current configuration, contruct neighbor list,
826 <       ! and calculate forces
767 <       call saveNeighborList(natoms, q)
825 >       !! save current configuration, construct neighbor list,
826 >       !! and calculate forces
827        
828 +       call saveNeighborList(nGroup, q_group)
829 +      
830         neighborListSize = size(list)
831 +       nlist = 0      
832        
833 <       nlist = 0
772 <      
773 <       do i = 1, natoms-1
833 >       do i = 1, nGroup-1
834            point(i) = nlist + 1
835            
836 <          inner: do j = i+1, natoms
836 >          do j = i+1, nGroup
837              
838 <             if (skipThisPair(i,j))  cycle inner
838 >             call get_interatomic_vector(q_group(:,i), &
839 >                  q_group(:,j), d_grp, rgrpsq)
840              
841 <             if (SIM_uses_molecular_cutoffs) then
781 <                call get_interatomic_vector(qcom(:,i), qcom(:,j), &
782 <                     dc, rcijsq)
783 <             else
784 <                call get_interatomic_vector(q(:,i), q(:,j), &
785 <                     dc, rcijsq)
786 <             endif
787 <            
788 <             if (rcijsq < rlistsq) then
789 <                
841 >             if (rgrpsq < rlistsq) then
842                  nlist = nlist + 1
843                  
844                  if (nlist > neighborListSize) then
845 <                   call expandNeighborList(natoms, listerror)
845 >                   call expandNeighborList(nGroup, listerror)
846                     if (listerror /= 0) then
847                        error = -1
848                        write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
# Line 800 | Line 852 | contains
852                  endif
853                  
854                  list(nlist) = j
855 +
856 +                vab = 0.0d0
857 +                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
858 +                     in_switching_region)
859                  
860 <                if (SIM_uses_molecular_cutoffs) then
861 <                   call get_interatomic_vector(q(:,i), q(:,j), &
862 <                        d, rijsq)
863 <                else
864 <                   d(1:3) = dc(1:3)
865 <                   rijsq = rcijsq
860 >                do ia = groupStart(i), groupStart(i+1)-1
861 >                   atom1 = groupList(ia)
862 >                  
863 >                   inner: do jb = groupStart(j), groupStart(j+1)-1
864 >                      atom2 = groupList(jb)
865 >                      
866 >                      if (skipThisPair(atom1, atom2)) cycle inner
867 >                      
868 >                      call get_interatomic_vector(q(:,atom1), &
869 >                           q(:,atom2), d_atm, ratmsq)
870 >                      
871 >                      call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
872 >                           do_pot, do_stress, &
873 >                           u_l, A, f, t, pot, vpair)
874 >                      
875 >                      vab = vab + vpair
876 >                      
877 >                   enddo inner
878 >                enddo
879 >
880 >                if (in_switching_region) then
881 >                   swderiv = vab*dswdr/rgrp
882 >                   do ia=groupStart(i), groupStart(i+1)-1
883 >                      atom1=groupList(ia)
884 >                      mf = mfact(atom1)                  
885 >                      f(1,atom1) = f(1,atom1) - swderiv*d_grp(1)*mf
886 >                      f(2,atom1) = f(2,atom1) - swderiv*d_grp(2)*mf
887 >                      f(3,atom1) = f(3,atom1) - swderiv*d_grp(3)*mf
888 >                   enddo
889 >                  
890 >                   do jb=groupStart(j), groupStart(j+1)-1
891 >                      atom2=groupList(jb)
892 >                      mf = mfact(atom2)
893 >                      f(1,atom2) = f(1,atom2) + swderiv*d_grp(1)*mf
894 >                      f(2,atom2) = f(2,atom2) + swderiv*d_grp(2)*mf
895 >                      f(3,atom2) = f(3,atom2) + swderiv*d_grp(3)*mf
896 >                   enddo
897                  endif
898 <                
899 <                call do_pair(i, j, rijsq, d, rcijsq, dc, mfact, &
900 <                     do_pot, do_stress, u_l, A, f, t, pot)
814 <                
815 <             endif
816 <          enddo inner
898 >
899 >             end if
900 >          enddo
901         enddo
902 +       point(nGroup) = nlist + 1          
903        
904 <       point(natoms) = nlist + 1
904 >    else  !! (of update_check)
905        
821    else !! (update)
822      
906         ! use the list to find the neighbors
907 <       do i = 1, natoms-1
907 >       do i = 1, nGroup-1
908            JBEG = POINT(i)
909            JEND = POINT(i+1) - 1
910 <          ! check thiat molecule i has neighbors
910 >          ! check that group i has neighbors
911            if (jbeg .le. jend) then
912              
913               do jnab = jbeg, jend
914                  j = list(jnab)
915  
916 +                call get_interatomic_vector(q_group(:,i), &
917 +                     q_group(:,j), d_grp, rgrpsq)
918  
919 <                call get_interatomic_vector(q(:,i), q(:,j), d, rijsq)
920 <                if (SIM_uses_molecular_cutoffs) then
921 <                   call get_interatomic_vector(qcom(:,i), qcom(:,j), &
922 <                        dc, rcijsq)
923 <                else
924 <                   dc(1:3) = d(1:3)
925 <                   rcijsq = rijsq
926 <                endif
927 <                
928 <                call do_pair(i, j, rijsq, d, rcijsq, dc, mfact, &
929 <                     do_pot, do_stress, u_l, A, f, t, pot)
919 >                vab = 0.0d0                    
920 >                call get_switch(rgrpsq, sw, dswdr, rgrp, group_switch, &
921 >                     in_switching_region)
922 >                                
923 >                do ia = groupStart(i), groupStart(i+1)-1
924 >                   atom1 = groupList(ia)
925 >                  
926 >                   do jb = groupStart(j), groupStart(j+1)-1
927 >                      atom2 = groupList(jb)                        
928 >                      
929 >                      call get_interatomic_vector(q(:,atom1), &
930 >                           q(:,atom2), d_atm, ratmsq)
931 >                      
932 >                      call do_pair(atom1, atom2, ratmsq, d_atm, sw, &
933 >                           do_pot, do_stress, &
934 >                           u_l, A, f, t, pot, vpair)
935 >                      
936 >                      vab = vab + vpair
937 >                      
938 >                   enddo
939 >                enddo
940                  
941 +                if (in_switching_region) then
942 +                   swderiv = vab*dswdr/rgrp
943 +                  
944 +                   do ia=groupStart(i), groupStart(i+1)-1
945 +                      atom1=groupList(ia)
946 +                      mf = mfact(atom1)                  
947 +                      f(1,atom1) = f(1,atom1) - swderiv*d_grp(1)*mf
948 +                      f(2,atom1) = f(2,atom1) - swderiv*d_grp(2)*mf
949 +                      f(3,atom1) = f(3,atom1) - swderiv*d_grp(3)*mf
950 +                   enddo
951 +                  
952 +                   do jb=groupStart(j), groupStart(j+1)-1
953 +                      atom2=groupList(jb)
954 +                      mf = mfact(atom2)
955 +                      f(1,atom2) = f(1,atom2) + swderiv*d_grp(1)*mf
956 +                      f(2,atom2) = f(2,atom2) + swderiv*d_grp(2)*mf
957 +                      f(3,atom2) = f(3,atom2) + swderiv*d_grp(3)*mf
958 +                   enddo
959 +                endif
960               enddo
961            endif
962         enddo
963      endif
964 <    
964 >        
965   #endif
966      
967      ! phew, done with main loop.
# Line 856 | Line 970 | contains
970   #ifdef PROFILE
971      call cpu_time(forceTimeFinal)
972      forceTime = forceTime + forceTimeFinal - forceTimeInitial
973 < #endif
973 > #endif    
974      
861    
975   #ifdef IS_MPI
976      !!distribute forces
977      
# Line 904 | Line 1017 | contains
1017         do i = 1, nlocal
1018            pot_local = pot_local + pot_Temp(i)
1019         enddo
1020 <
1020 >      
1021      endif
1022   #endif
1023      
# Line 977 | Line 1090 | contains
1090      
1091      
1092    end subroutine do_force_loop
1093 +
1094    
1095 <  subroutine do_pair(i, j, rijsq, d, rcijsq, dc, mfact, do_pot, do_stress, &
1096 <       u_l, A, f, t, pot)
1095 >  subroutine do_pair(i, j, rijsq, d, sw, do_pot, do_stress, &
1096 >       u_l, A, f, t, pot, vpair)
1097  
1098 <    real( kind = dp ) :: pot
1098 >    real( kind = dp ) :: pot, vpair, sw
1099      real( kind = dp ), dimension(nLocal)   :: mfact
1100      real( kind = dp ), dimension(3,nLocal) :: u_l
1101      real( kind = dp ), dimension(9,nLocal) :: A
# Line 990 | Line 1104 | contains
1104  
1105      logical, intent(inout) :: do_pot, do_stress
1106      integer, intent(in) :: i, j
1107 <    real ( kind = dp ), intent(inout) :: rijsq, rcijsq
1108 <    real ( kind = dp )                :: r, rc
1109 <    real ( kind = dp ), intent(inout) :: d(3), dc(3)
1107 >    real ( kind = dp ), intent(inout) :: rijsq
1108 >    real ( kind = dp )                :: r
1109 >    real ( kind = dp ), intent(inout) :: d(3)
1110      integer :: me_i, me_j
1111  
1112      r = sqrt(rijsq)
999    if (SIM_uses_molecular_cutoffs) then
1000       rc = sqrt(rcijsq)
1001    else
1002       rc = r
1003    endif
1113  
1005
1114   #ifdef IS_MPI
1115      if (tagRow(i) .eq. tagColumn(j)) then
1116         write(0,*) 'do_pair is doing', i , j, tagRow(i), tagColumn(j)
# Line 1017 | Line 1125 | contains
1125      if (FF_uses_LJ .and. SIM_uses_LJ) then
1126        
1127         if ( PropertyMap(me_i)%is_LJ .and. PropertyMap(me_j)%is_LJ ) then
1128 <          call do_lj_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
1128 >          call do_lj_pair(i, j, d, r, rijsq, sw, vpair, pot, f, do_pot, &
1129 >               do_stress)
1130         endif
1131        
1132      endif
# Line 1025 | Line 1134 | contains
1134      if (FF_uses_charges .and. SIM_uses_charges) then
1135        
1136         if (PropertyMap(me_i)%is_Charge .and. PropertyMap(me_j)%is_Charge) then
1137 <          call do_charge_pair(i, j, d, r, rijsq, dc, rc, rcijsq, mfact, &
1138 <               pot, f, do_pot, do_stress, SIM_uses_molecular_cutoffs)
1137 >          call do_charge_pair(i, j, d, r, rijsq, sw, vpair, pot, f, do_pot, &
1138 >               do_stress)
1139         endif
1140        
1141      endif
# Line 1034 | Line 1143 | contains
1143      if (FF_uses_dipoles .and. SIM_uses_dipoles) then
1144        
1145         if ( PropertyMap(me_i)%is_DP .and. PropertyMap(me_j)%is_DP) then
1146 <          call do_dipole_pair(i, j, d, r, rijsq, pot, u_l, f, t, &
1146 >          call do_dipole_pair(i, j, d, r, rijsq, sw, vpair, pot, u_l, f, t, &
1147                 do_pot, do_stress)
1148            if (FF_uses_RF .and. SIM_uses_RF) then
1149               call accumulate_rf(i, j, r, u_l)
# Line 1047 | Line 1156 | contains
1156      if (FF_uses_Sticky .and. SIM_uses_sticky) then
1157  
1158         if ( PropertyMap(me_i)%is_Sticky .and. PropertyMap(me_j)%is_Sticky) then
1159 <          call do_sticky_pair(i, j, d, r, rijsq, A, pot, f, t, &
1159 >          call do_sticky_pair(i, j, d, r, rijsq, sw, vpair, pot, A, f, t, &
1160                 do_pot, do_stress)
1161         endif
1162  
# Line 1057 | Line 1166 | contains
1166      if (FF_uses_GB .and. SIM_uses_GB) then
1167        
1168         if ( PropertyMap(me_i)%is_GB .and. PropertyMap(me_j)%is_GB) then
1169 <          call do_gb_pair(i, j, d, r, rijsq, u_l, pot, f, t, &
1169 >          call do_gb_pair(i, j, d, r, rijsq, sw, vpair, pot, u_l, f, t, &
1170                 do_pot, do_stress)          
1171         endif
1172  
# Line 1066 | Line 1175 | contains
1175      if (FF_uses_EAM .and. SIM_uses_EAM) then
1176        
1177         if ( PropertyMap(me_i)%is_EAM .and. PropertyMap(me_j)%is_EAM) then
1178 <          call do_eam_pair(i, j, d, r, rijsq, pot, f, do_pot, do_stress)
1178 >          call do_eam_pair(i, j, d, r, rijsq, sw, vpair, pot, f, &
1179 >               do_pot, do_stress)
1180         endif
1181        
1182      endif
1183 <
1183 >    
1184    end subroutine do_pair
1075
1076
1185  
1186    subroutine do_prepair(i, j, rijsq, d, rcijsq, dc, &
1187         do_pot, do_stress, u_l, A, f, t, pot)
# Line 1199 | Line 1307 | contains
1307     q_Row = 0.0_dp
1308     q_Col = 0.0_dp
1309  
1310 <   qcom_Row = 0.0_dp
1311 <   qcom_Col = 0.0_dp  
1310 >   q_group_Row = 0.0_dp
1311 >   q_group_Col = 0.0_dp  
1312    
1313     u_l_Row = 0.0_dp
1314     u_l_Col = 0.0_dp

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines