ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE_old/src/mdtools/libmdCode/do_Forces.F90
(Generate patch)

Comparing trunk/OOPSE_old/src/mdtools/libmdCode/do_Forces.F90 (file contents):
Revision 293 by mmeineke, Thu Mar 6 16:07:57 2003 UTC vs.
Revision 297 by gezelter, Thu Mar 6 22:08:29 2003 UTC

# Line 1 | Line 1
1   !! Calculates Long Range forces.
2   !! @author Charles F. Vardeman II
3   !! @author Matthew Meineke
4 < !! @version $Id: do_Forces.F90,v 1.2 2003-03-06 16:07:55 mmeineke Exp $, $Date: 2003-03-06 16:07:55 $, $Name: not supported by cvs2svn $, $Revision: 1.2 $
4 > !! @version $Id: do_Forces.F90,v 1.4 2003-03-06 22:08:29 gezelter Exp $, $Date: 2003-03-06 22:08:29 $, $Name: not supported by cvs2svn $, $Revision: 1.4 $
5  
6  
7  
# Line 30 | Line 30 | module do_Forces
30    type (atype), pointer :: ListTail => null()
31  
32  
33 – !! LJ mixing array  
34 – !  type (lj_atype), dimension(:,:), pointer :: ljMixed => null()
33  
34  
35    logical, save :: firstTime = .True.
# Line 50 | Line 48 | module do_Forces
48   !! Logical has lj force field module been initialized?
49    logical, save :: isFFinit = .false.
50  
51 + !! Use periodic boundry conditions
52 +  logical :: wrap = .false.
53  
54 + !! Potential energy global module variables
55 + #ifdef IS_MPI
56 +  real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow = 0.0_dp
57 +  real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol = 0.0_dp
58 +
59 +  real(kind = dp), dimension(3,getNrow(plan_row)) :: muRow = 0.0_dp
60 +  real(kind = dp), dimension(3,getNcol(plan_col)) :: muCol = 0.0_dp
61 +
62 +  real(kind = dp), dimension(3,getNrow(plan_row)) :: u_lRow = 0.0_dp
63 +  real(kind = dp), dimension(3,getNcol(plan_col)) :: u_lCol = 0.0_dp
64 +
65 +  real(kind = dp), dimension(9,getNrow(plan_row)) :: ARow = 0.0_dp
66 +  real(kind = dp), dimension(9,getNcol(plan_col)) :: ACol = 0.0_dp  
67 +
68 +  real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow = 0.0_dp
69 +  real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol = 0.0_dp
70 +  real(kind = dp), dimension(3,getNlocal()) :: fTemp1 = 0.0_dp
71 +  real(kind = dp), dimension(3,getNlocal()) :: tTemp1 = 0.0_dp
72 +  real(kind = dp), dimension(3,getNlocal()) :: fTemp2 = 0.0_dp
73 +  real(kind = dp), dimension(3,getNlocal()) :: tTemp2 = 0.0_dp
74 +  real(kind = dp), dimension(3,getNlocal()) :: fTemp = 0.0_dp
75 +  real(kind = dp), dimension(3,getNlocal()) :: tTemp = 0.0_dp
76 +
77 +  real(kind = dp), dimension(3,getNrow(plan_row)) :: tRow = 0.0_dp
78 +  real(kind = dp), dimension(3,getNcol(plan_col)) :: tCol = 0.0_dp
79 +
80 +  real(kind = dp), dimension(3,getNrow(plan_row)) :: rflRow = 0.0_dp
81 +  real(kind = dp), dimension(3,getNcol(plan_col)) :: rflCol = 0.0_dp
82 +  real(kind = dp), dimension(3,getNlocal()) :: rflTemp = 0.0_dp
83 +
84 +  real(kind = dp), dimension(getNrow(plan_row)) :: eRow = 0.0_dp
85 +  real(kind = dp), dimension(getNcol(plan_col)) :: eCol = 0.0_dp
86 +
87 +  real(kind = dp), dimension(getNlocal()) :: eTemp = 0.0_dp
88 + #endif
89 +  real(kind = dp) :: pe = 0.0_dp
90 +  real(kind = dp), dimension(3,natoms) :: fTemp = 0.0_dp
91 +  real(kind = dp), dimension(3,natoms) :: tTemp = 0.0_dp
92 +  real(kind = dp), dimension(3,natoms) :: rflTemp = 0.0_dp
93 +  real(kind = dp), dimension(9) :: tauTemp = 0.0_dp
94 +
95 +  logical :: do_preForce  = .false.
96 +  logical :: do_postForce = .false.
97 +
98 +
99 +
100   !! Public methods and data
101    public :: new_atype
102    public :: do_forceLoop
57 –  public :: getLjPot
103    public :: init_FF
104  
105    
# Line 258 | Line 303 | contains
303  
304   !! FORCE routine Calculates Lennard Jones forces.
305   !------------------------------------------------------------->
306 <  subroutine do__force_loop(q,A,mu,u_l,f,t,tau,potE,do_pot,FFerror)
306 >  subroutine do_force_loop(q,A,mu,u_l,f,t,tau,potE,do_pot,FFerror)
307   !! Position array provided by C, dimensioned by getNlocal
308      real ( kind = dp ), dimension(3,getNlocal()) :: q
309    !! Rotation Matrix for each long range particle in simulation.
# Line 275 | Line 320 | contains
320  
321   !! Stress Tensor
322      real( kind = dp), dimension(9) :: tau
323 <    real( kind = dp), dimension(9) :: tauTemp
323 >
324      real ( kind = dp ) :: potE
325      logical ( kind = 2) :: do_pot
326      integer :: FFerror
# Line 293 | Line 338 | contains
338    real( kind = DP ) :: pot_local
339  
340   !! Local arrays needed for MPI
296 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow = 0.0_dp
297 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol = 0.0_dp
341  
299 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: muRow = 0.0_dp
300 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: muCol = 0.0_dp
301 –
302 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: u_lRow = 0.0_dp
303 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: u_lCol = 0.0_dp
304 –
305 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: ARow = 0.0_dp
306 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: ACol = 0.0_dp
307 –
308 –  
309 –
310 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow = 0.0_dp
311 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol = 0.0_dp
312 –  real(kind = dp), dimension(3,getNlocal()) :: fMPITemp = 0.0_dp
313 –
314 –  real(kind = dp), dimension(3,getNrow(plan_row)) :: tRow = 0.0_dp
315 –  real(kind = dp), dimension(3,getNcol(plan_col)) :: tCol = 0.0_dp
316 –  real(kind = dp), dimension(3,getNlocal()) :: tMPITemp = 0.0_dp
317 –
318 –
319 –  real(kind = dp), dimension(getNrow(plan_row)) :: eRow = 0.0_dp
320 –  real(kind = dp), dimension(getNcol(plan_col)) :: eCol = 0.0_dp
321 –
322 –  real(kind = dp), dimension(getNlocal()) :: eTemp = 0.0_dp
323 –
342   #endif
343  
344  
# Line 332 | Line 350 | contains
350    integer ::  i, j, jbeg, jend, jnab, idim, jdim, idim2, jdim2, dim, dim2
351    integer :: nlist
352    integer :: j_start
353 <  integer :: tag_i,tag_j
354 <  real( kind = DP ) ::  r, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp
355 <  real( kind = dp ) :: fx,fy,fz
356 <  real( kind = DP ) ::  drdx, drdy, drdz
339 <  real( kind = DP ) ::  rxi, ryi, rzi, rxij, ryij, rzij, rijsq
353 >
354 >  real( kind = DP ) ::  r_ij, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp
355 >
356 >  real( kind = DP ) ::  rx_ij, ry_ij, rz_ij, rijsq
357    real( kind = DP ) ::  rlistsq, rcutsq,rlist,rcut
358  
359    real( kind = DP ) :: dielectric = 0.0_dp
360  
361   ! a rig that need to be fixed.
362   #ifdef IS_MPI
346 –  logical :: newtons_thrd = .true.
363    real( kind = dp ) :: pe_local
364    integer :: nlocal
365   #endif
# Line 376 | Line 392 | contains
392    natoms = getNlocal()
393    call getRcut(rcut,rcut2=rcutsq)
394    call getRlist(rlist,rlistsq)
395 +
396   !! Find ensemble
397    if (isEnsemble("NPT")) do_stress = .true.
398 + !! set to wrap
399 +  if (isPBC()) wrap = .true.
400  
401 < #ifndef IS_MPI
383 <  nrow = natoms - 1
384 <  ncol = natoms
385 < #else
386 <  nrow = getNrow(plan_row)
387 <  ncol = getNcol(plan_col)
388 <  nlocal = natoms
389 <  j_start = 1
390 < #endif
401 >
402  
403    
404   !! See if we need to update neighbor lists
405    call check(q,update_nlist)
395 – !  if (firstTime) then
396 – !     update_nlist = .true.
397 – !     firstTime = .false.
398 – !  endif
406  
407   !--------------WARNING...........................
408   ! Zero variables, NOTE:::: Forces are zeroed in C
409   ! Zeroing them here could delete previously computed
410   ! Forces.
411   !------------------------------------------------
412 < #ifndef IS_MPI
406 < !  nloops = nloops + 1
407 <  pe = 0.0E0_DP
408 <
409 < #else
410 <    fRow = 0.0E0_DP
411 <    fCol = 0.0E0_DP
412 >  call zero_module_variables()
413  
413 –    pe_local = 0.0E0_DP
414  
415 –    eRow = 0.0E0_DP
416 –    eCol = 0.0E0_DP
417 –    eTemp = 0.0E0_DP
418 – #endif
419 –
415   ! communicate MPI positions
416   #ifdef IS_MPI    
417      call gather(q,qRow,plan_row3d)
# Line 430 | Line 425 | contains
425  
426      call gather(A,ARow,plan_row_rotation)
427      call gather(A,ACol,plan_col_rotation)
433 –
428   #endif
429  
430  
431 <  if (update_nlist) then
438 <
439 <     ! save current configuration, contruct neighbor list,
440 <     ! and calculate forces
441 <     call save_neighborList(q)
442 <    
443 <     neighborListSize = getNeighborListSize()
444 <     nlist = 0
445 <    
431 > #ifdef IS_MPI
432      
433 +    if (update_nlist) then
434 +      
435 +       ! save current configuration, contruct neighbor list,
436 +       ! and calculate forces
437 +       call save_neighborList(q)
438 +      
439 +       neighborListSize = getNeighborListSize()
440 +       nlist = 0
441 +      
442 +       nrow = getNrow(plan_row)
443 +       ncol = getNcol(plan_col)
444 +       nlocal = getNlocal()
445 +      
446 +       do i = 1, nrow
447 +          point(i) = nlist + 1
448 +          Atype_i => identPtrListRow(i)%this
449 +          
450 +          inner: do j = 1, ncol
451 +             Atype_j => identPtrListColumn(j)%this
452 +            
453 +             call get_interatomic_vector(i,j,qRow(:,i),qCol(:,j),&
454 +                  rxij,ryij,rzij,rijsq,r)
455 +            
456 +             ! skip the loop if the atoms are identical
457 +             if (mpi_cycle_jLoop(i,j)) cycle inner:
458 +            
459 +             if (rijsq <  rlistsq) then            
460 +                
461 +                nlist = nlist + 1
462 +                
463 +                if (nlist > neighborListSize) then
464 +                   call expandList(listerror)
465 +                   if (listerror /= 0) then
466 +                      FFerror = -1
467 +                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
468 +                      return
469 +                   end if
470 +                endif
471 +                
472 +                list(nlist) = j
473 +                
474 +                
475 +                if (rijsq <  rcutsq) then
476 +                   call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
477 +                endif
478 +             endif
479 +          enddo inner
480 +       enddo
481  
482 <     do i = 1, nrow
483 <        point(i) = nlist + 1
484 < #ifdef IS_MPI
485 <        ljAtype_i => identPtrListRow(i)%this
486 <        tag_i = tagRow(i)
487 <        rxi = qRow(1,i)
488 <        ryi = qRow(2,i)
489 <        rzi = qRow(3,i)
482 >       point(nrow + 1) = nlist + 1
483 >      
484 >    else !! (update)
485 >
486 >       ! use the list to find the neighbors
487 >       do i = 1, nrow
488 >          JBEG = POINT(i)
489 >          JEND = POINT(i+1) - 1
490 >          ! check thiat molecule i has neighbors
491 >          if (jbeg .le. jend) then
492 >
493 >             Atype_i => identPtrListRow(i)%this
494 >             do jnab = jbeg, jend
495 >                j = list(jnab)
496 >                Atype_j = identPtrListColumn(j)%this
497 >                call get_interatomic_vector(i,j,qRow(:,i),qCol(:,j),&
498 >                     rxij,ryij,rzij,rijsq,r)
499 >                
500 >                call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
501 >             enddo
502 >          endif
503 >       enddo
504 >    endif
505 >
506   #else
507 <        ljAtype_i   => identPtrList(i)%this
508 <        j_start = i + 1
509 <        rxi = q(1,i)
510 <        ryi = q(2,i)
511 <        rzi = q(3,i)
507 >    
508 >    if (update_nlist) then
509 >      
510 >       ! save current configuration, contruct neighbor list,
511 >       ! and calculate forces
512 >       call save_neighborList(q)
513 >      
514 >       neighborListSize = getNeighborListSize()
515 >       nlist = 0
516 >      
517 >    
518 >       do i = 1, natoms-1
519 >          point(i) = nlist + 1
520 >          Atype_i   => identPtrList(i)%this
521 >
522 >          inner: do j = i+1, natoms
523 >             Atype_j   => identPtrList(j)%this
524 >             call get_interatomic_vector(i,j,q(:,i),q(:,j),&
525 >                  rxij,ryij,rzij,rijsq,r)
526 >          
527 >             if (rijsq <  rlistsq) then
528 >
529 >                nlist = nlist + 1
530 >                
531 >                if (nlist > neighborListSize) then
532 >                   call expandList(listerror)
533 >                   if (listerror /= 0) then
534 >                      FFerror = -1
535 >                      write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
536 >                      return
537 >                   end if
538 >                endif
539 >                
540 >                list(nlist) = j
541 >
542 >    
543 >                if (rijsq <  rcutsq) then
544 >                   call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
545 >                endif
546 >             endif
547 >          enddo inner
548 >       enddo
549 >      
550 >       point(natoms) = nlist + 1
551 >      
552 >    else !! (update)
553 >      
554 >       ! use the list to find the neighbors
555 >       do i = 1, nrow
556 >          JBEG = POINT(i)
557 >          JEND = POINT(i+1) - 1
558 >          ! check thiat molecule i has neighbors
559 >          if (jbeg .le. jend) then
560 >            
561 >             Atype_i => identPtrList(i)%this
562 >             do jnab = jbeg, jend
563 >                j = list(jnab)
564 >                Atype_j = identPtrList(j)%this
565 >                call get_interatomic_vector(i,j,q(:,i),q(:,j),&
566 >                     rxij,ryij,rzij,rijsq,r)
567 >                call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij)
568 >             enddo
569 >          endif
570 >       enddo
571 >    endif
572 >
573 > #endif
574 >
575 >
576 > #ifdef IS_MPI
577 >    !! distribute all reaction field stuff (or anything for post-pair):
578 >    call scatter(rflRow,rflTemp1,plan_row3d)
579 >    call scatter(rflCol,rflTemp2,plan_col3d)
580 >    do i = 1,nlocal
581 >       rflTemp(1:3,i) = rflTemp1(1:3,i) + rflTemp2(1:3,i)
582 >    end do
583   #endif
584  
585 <        inner: do j = j_start, ncol
585 > ! This is the post-pair loop:
586   #ifdef IS_MPI
466 – ! Assign identity pointers and tags
467 –           ljAtype_j => identPtrListColumn(j)%this
468 –           tag_j = tagColumn(j)
469 –           if (newtons_thrd) then
470 –              if (tag_i <= tag_j) then
471 –                 if (mod(tag_i + tag_j,2) == 0) cycle inner
472 –              else                
473 –                 if (mod(tag_i + tag_j,2) == 1) cycle inner
474 –              endif
475 –           endif
587  
588 <           rxij = wrap(rxi - qCol(1,j), 1)
589 <           ryij = wrap(ryi - qCol(2,j), 2)
590 <           rzij = wrap(rzi - qCol(3,j), 3)
591 < #else          
592 <           ljAtype_j   => identPtrList(j)%this
593 <           rxij = wrap(rxi - q(1,j), 1)
483 <           ryij = wrap(ryi - q(2,j), 2)
484 <           rzij = wrap(rzi - q(3,j), 3)
485 <      
486 < #endif          
487 <           rijsq = rxij*rxij + ryij*ryij + rzij*rzij
588 >    if (system_has_postpair_atoms) then
589 >       do i = 1, nlocal
590 >          Atype_i => identPtrListRow(i)%this
591 >          call do_postpair(i, Atype_i)
592 >       enddo
593 >    endif
594  
595 + #else
596 +
597 +    if (system_has_postpair_atoms) then
598 +       do i = 1, natoms
599 +          Atype_i => identPtr(i)%this
600 +          call do_postpair(i, Atype_i)
601 +       enddo
602 +    endif
603 +
604 + #endif
605 +
606 +
607 +
608 +
609   #ifdef IS_MPI
610 <             if (rijsq <=  rlistsq .AND. &
611 <                  tag_j /= tag_i) then
610 >    !!distribute forces
611 >
612 >    call scatter(fRow,fTemp1,plan_row3d)
613 >    call scatter(fCol,fTemp2,plan_col3d)
614 >
615 >
616 >    do i = 1,nlocal
617 >       fTemp(1:3,i) = fTemp1(1:3,i) + fTemp2(1:3,i)
618 >    end do
619 >
620 >    if (do_torque) then
621 >       call scatter(tRow,tTemp1,plan_row3d)
622 >       call scatter(tCol,tTemp2,plan_col3d)
623 >    
624 >       do i = 1,nlocal
625 >          tTemp(1:3,i) = tTemp1(1:3,i) + tTemp2(1:3,i)
626 >       end do
627 >    endif
628 >
629 >    if (do_pot) then
630 >       ! scatter/gather pot_row into the members of my column
631 >       call scatter(eRow,eTemp,plan_row)
632 >      
633 >       ! scatter/gather pot_local into all other procs
634 >       ! add resultant to get total pot
635 >       do i = 1, nlocal
636 >          pe_local = pe_local + eTemp(i)
637 >       enddo
638 >
639 >       eTemp = 0.0E0_DP
640 >       call scatter(eCol,eTemp,plan_col)
641 >       do i = 1, nlocal
642 >          pe_local = pe_local + eTemp(i)
643 >       enddo
644 >      
645 >       pe = pe_local
646 >    endif
647   #else
648 <          
649 <             if (rijsq <  rlistsq) then
648 > ! Copy local array into return array for c
649 >    f = f+fTemp
650 >    t = t+tTemp
651   #endif
496 –            
497 –              nlist = nlist + 1
498 –              if (nlist > neighborListSize) then
499 –                 call expandList(listerror)
500 –                 if (listerror /= 0) then
501 –                    FFerror = -1
502 –                    write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded."
503 –                    return
504 –                 end if
505 –              endif
506 –              list(nlist) = j
652  
653 +    potE = pe
654 +
655 +
656 +    if (do_stress) then
657 + #ifdef IS_MPI
658 +       mpi_allreduce = (tau,tauTemp,9,mpi_double_precision,mpi_sum, &
659 +            mpi_comm_world,mpi_err)
660 + #else
661 +       tau = tauTemp
662 + #endif      
663 +    endif
664 +
665 +  end subroutine do_force_loop
666 +
667 +
668 +
669 +
670 +
671 +
672 +
673 +
674 +
675 +
676 + !! Calculate any pre-force loop components and update nlist if necessary.
677 +  subroutine do_preForce(updateNlist)
678 +    logical, intent(inout) :: updateNlist
679 +
680 +
681 +
682 +  end subroutine do_preForce
683 +
684 +
685 +
686 +
687 +
688 +
689 +
690 +
691 +
692 +
693 +
694 +
695 +
696 + !! Calculate any post force loop components, i.e. reaction field, etc.
697 +  subroutine do_postForce()
698 +
699 +
700 +
701 +  end subroutine do_postForce
702 +
703 +
704 +
705 +
706 +
707 +
708 +
709 +
710 +
711 +
712 +
713 +
714 +
715 +
716 +
717 +
718 +
719 +  subroutine do_pair(atype_i,atype_j,i,j,r_ij,rx_ij,ry_ij,rz_ij)
720 +    type (atype ), pointer, intent(inout) :: atype_i
721 +    type (atype ), pointer, intent(inout) :: atype_j
722 +    integer :: i
723 +    integer :: j
724 +    real ( kind = dp ), intent(inout) :: rx_ij
725 +    real ( kind = dp ), intent(inout) :: ry_ij
726 +    real ( kind = dp ), intent(inout) :: rz_ij
727 +
728 +
729 +    real( kind = dp ) :: fx = 0.0_dp
730 +    real( kind = dp ) :: fy = 0.0_dp
731 +    real( kind = dp ) :: fz = 0.0_dp  
732 +  
733 +    real( kind = dp ) ::  drdx = 0.0_dp
734 +    real( kind = dp ) ::  drdy = 0.0_dp
735 +    real( kind = dp ) ::  drdz = 0.0_dp
736      
509 –              if (rijsq <  rcutsq) then
510 –                
511 –                 r = dsqrt(rijsq)
512 –      
513 –                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
514 –      
515 – #ifdef IS_MPI
516 –                eRow(i) = eRow(i) + pot*0.5
517 –                eCol(i) = eCol(i) + pot*0.5
518 – #else
519 –                    pe = pe + pot
520 – #endif                
521 –            
522 –                drdx = -rxij / r
523 –                drdy = -ryij / r
524 –                drdz = -rzij / r
525 –                
526 –                fx = dudr * drdx
527 –                fy = dudr * drdy
528 –                fz = dudr * drdz
529 –                
530 – #ifdef IS_MPI
531 –                fCol(1,j) = fCol(1,j) - fx
532 –                fCol(2,j) = fCol(2,j) - fy
533 –                fCol(3,j) = fCol(3,j) - fz
534 –                
535 –                fRow(1,j) = fRow(1,j) + fx
536 –                fRow(2,j) = fRow(2,j) + fy
537 –                fRow(3,j) = fRow(3,j) + fz
538 – #else
539 –                f(1,j) = f(1,j) - fx
540 –                f(2,j) = f(2,j) - fy
541 –                f(3,j) = f(3,j) - fz
542 –                f(1,i) = f(1,i) + fx
543 –                f(2,i) = f(2,i) + fy
544 –                f(3,i) = f(3,i) + fz
545 – #endif
546 –                
547 –                if (do_stress) then
548 –                   tauTemp(1) = tauTemp(1) + fx * rxij
549 –                   tauTemp(2) = tauTemp(2) + fx * ryij
550 –                   tauTemp(3) = tauTemp(3) + fx * rzij
551 –                   tauTemp(4) = tauTemp(4) + fy * rxij
552 –                   tauTemp(5) = tauTemp(5) + fy * ryij
553 –                   tauTemp(6) = tauTemp(6) + fy * rzij
554 –                   tauTemp(7) = tauTemp(7) + fz * rxij
555 –                   tauTemp(8) = tauTemp(8) + fz * ryij
556 –                   tauTemp(9) = tauTemp(9) + fz * rzij
557 –                endif
558 –             endif
559 –          enddo inner
560 –     enddo
737  
738 < #ifdef IS_MPI
739 <     point(nrow + 1) = nlist + 1
740 < #else
565 <     point(natoms) = nlist + 1
566 < #endif
738 >    if (Atype_i%is_LJ .and. Atype_j%is_LJ) then
739 >       call getLJForce(r,pot,dudr,ljAtype_i,ljAtype_j,fx,fy,fz)
740 >    endif
741  
742 <  else
742 >    if (Atype_i%is_dp .and. Atype_j%is_dp) then
743  
570 –     ! use the list to find the neighbors
571 –     do i = 1, nrow
572 –        JBEG = POINT(i)
573 –        JEND = POINT(i+1) - 1
574 –        ! check thiat molecule i has neighbors
575 –        if (jbeg .le. jend) then
744   #ifdef IS_MPI
745 <           ljAtype_i => identPtrListRow(i)%this
746 <           rxi = qRow(1,i)
579 <           ryi = qRow(2,i)
580 <           rzi = qRow(3,i)
745 >       call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
746 >            ulRow(:,i), ulCol(:,j), rt, rrf, pot)
747   #else
748 <           ljAtype_i   => identPtrList(i)%this
749 <           rxi = q(1,i)
584 <           ryi = q(2,i)
585 <           rzi = q(3,i)
748 >       call dipole_dipole(i, j, atype_i, atype_j, rx_ij, ry_ij, rz_ij, r_ij, &
749 >            ul(:,i), ul(:,j), rt, rrf, pot)
750   #endif
751 <           do jnab = jbeg, jend
752 <              j = list(jnab)
751 >
752 >       if (do_reaction_field) then
753   #ifdef IS_MPI
754 <              ljAtype_j = identPtrListColumn(j)%this
755 <              rxij = wrap(rxi - qCol(1,j), 1)
592 <              ryij = wrap(ryi - qCol(2,j), 2)
593 <              rzij = wrap(rzi - qCol(3,j), 3)
754 >          call accumulate_rf(i, j, r_ij, rflRow(:,i), rflCol(:j), &
755 >               ulRow(:i), ulCol(:,j), rt, rrf)
756   #else
757 <              ljAtype_j = identPtrList(j)%this
758 <              rxij = wrap(rxi - q(1,j), 1)
597 <              ryij = wrap(ryi - q(2,j), 2)
598 <              rzij = wrap(rzi - q(3,j), 3)
757 >          call accumulate_rf(i, j, r_ij, rfl(:,i), rfl(:j), &
758 >               ul(:,i), ul(:,j), rt, rrf)
759   #endif
760 <              rijsq = rxij*rxij + ryij*ryij + rzij*rzij
601 <              
602 <              if (rijsq <  rcutsq) then
760 >       endif
761  
762 <                 r = dsqrt(rijsq)
763 <                
764 <                 call getLJPot(r,pot,dudr,ljAtype_i,ljAtype_j)
762 >
763 >    endif
764 >
765 >    if (Atype_i%is_sticky .and. Atype_j%is_sticky) then
766 >       call getstickyforce(r,pot,dudr,ljAtype_i,ljAtype_j)
767 >    endif
768 >
769 >      
770   #ifdef IS_MPI
771                  eRow(i) = eRow(i) + pot*0.5
772                  eCol(i) = eCol(i) + pot*0.5
773   #else
774 <                pe = pe + pot
774 >                    pe = pe + pot
775   #endif                
776 <  
776 >            
777                  drdx = -rxij / r
778                  drdy = -ryij / r
779                  drdz = -rzij / r
# Line 618 | Line 781 | contains
781                  fx = dudr * drdx
782                  fy = dudr * drdy
783                  fz = dudr * drdz
784 +
785 +
786 +
787 +
788 +
789 +
790                  
791   #ifdef IS_MPI
792                  fCol(1,j) = fCol(1,j) - fx
# Line 628 | Line 797 | contains
797                  fRow(2,j) = fRow(2,j) + fy
798                  fRow(3,j) = fRow(3,j) + fz
799   #else
800 <                f(1,j) = f(1,j) - fx
801 <                f(2,j) = f(2,j) - fy
802 <                f(3,j) = f(3,j) - fz
803 <                f(1,i) = f(1,i) + fx
804 <                f(2,i) = f(2,i) + fy
805 <                f(3,i) = f(3,i) + fz
800 >                fTemp(1,j) = fTemp(1,j) - fx
801 >                fTemp(2,j) = fTemp(2,j) - fy
802 >                fTemp(3,j) = fTemp(3,j) - fz
803 >                fTemp(1,i) = fTemp(1,i) + fx
804 >                fTemp(2,i) = fTemp(2,i) + fy
805 >                fTemp(3,i) = fTemp(3,i) + fz
806   #endif
807                  
808                  if (do_stress) then
# Line 647 | Line 816 | contains
816                     tauTemp(8) = tauTemp(8) + fz * ryij
817                     tauTemp(9) = tauTemp(9) + fz * rzij
818                  endif
650 –                
651 –                
652 –             endif
653 –          enddo
654 –       endif
655 –    enddo
656 – endif
657 –
819  
820  
660 – #ifdef IS_MPI
661 –    !!distribute forces
821  
822 <    call scatter(fRow,f,plan_row3d)
664 <    call scatter(fCol,fMPITemp,plan_col3d)
822 >  end subroutine do_pair
823  
666 –    do i = 1,nlocal
667 –       f(1:3,i) = f(1:3,i) + fMPITemp(1:3,i)
668 –    end do
824  
825  
826 <    call scatter(tRow,t,plan_row3d)
827 <    call scatter(tCol,tMPITemp,plan_col3d)
826 >
827 >
828 >
829 >
830 >
831 >
832 >
833 >
834 >
835 >
836 >
837 >
838 >
839 >  subroutine get_interatomic_vector(q_i,q_j,rx_ij,ry_ij,rz_ij,r_sq,r_ij)
840 > !---------------- Arguments-------------------------------
841 >   !! index i
842 >
843 >    !! Position array
844 >    real (kind = dp), dimension(3) :: q_i
845 >    real (kind = dp), dimension(3) :: q_j
846 >    !! x component of vector between i and j
847 >    real ( kind = dp ), intent(out)  :: rx_ij
848 >    !! y component of vector between i and j
849 >    real ( kind = dp ), intent(out)  :: ry_ij
850 >    !! z component of vector between i and j
851 >    real ( kind = dp ), intent(out)  :: rz_ij
852 >    !! magnitude of r squared
853 >    real ( kind = dp ), intent(out) :: r_sq
854 >    !! magnitude of vector r between atoms i and j.
855 >    real ( kind = dp ), intent(out) :: r_ij
856 >    !! wrap into periodic box.
857 >    logical, intent(in) :: wrap
858 >
859 > !--------------- Local Variables---------------------------
860 >    !! Distance between i and j
861 >    real( kind = dp ) :: d(3)
862 > !---------------- END DECLARATIONS-------------------------
863 >
864 >
865 > ! Find distance between i and j
866 >    d(1:3) = q_i(1:3) - q_j(1:3)
867 >
868 > ! Wrap back into periodic box if necessary
869 >    if ( wrap ) then
870 >       d(1:3) = d(1:3) - thisSim%box(1:3) * sign(1.0_dp,thisSim%box(1:3)) * &
871 >            int(abs(d(1:3)/thisSim%box(1:3) + 0.5_dp)
872 >    end if
873      
874 <    do i = 1,nlocal
875 <       t(1:3,i) = t(1:3,i) + tMPITemp(1:3,i)
876 <    end do
874 > !   Find Magnitude of the vector
875 >    r_sq = dot_product(d,d)
876 >    r_ij = sqrt(r_sq)
877  
878 + !   Set each component for force calculation
879 +    rx_ij = d(1)
880 +    ry_ij = d(2)
881 +    rz_ij = d(3)
882  
679 –    if (do_pot) then
680 –       ! scatter/gather pot_row into the members of my column
681 –       call scatter(eRow,eTemp,plan_row)
682 –      
683 –       ! scatter/gather pot_local into all other procs
684 –       ! add resultant to get total pot
685 –       do i = 1, nlocal
686 –          pe_local = pe_local + eTemp(i)
687 –       enddo
688 –       if (newtons_thrd) then
689 –          eTemp = 0.0E0_DP
690 –          call scatter(eCol,eTemp,plan_col)
691 –          do i = 1, nlocal
692 –             pe_local = pe_local + eTemp(i)
693 –          enddo
694 –       endif
695 –       pe = pe_local
696 –    endif
883  
884 < #endif
884 >  end subroutine get_interatomic_vector
885  
886 <    potE = pe
886 >  subroutine zero_module_variables()
887  
888 + #ifndef IS_MPI
889  
890 <    if (do_stress) then
891 < #ifdef IS_MPI
892 <       mpi_allreduce = (tau,tauTemp,9,mpi_double_precision,mpi_sum, &
893 <            mpi_comm_world,mpi_err)
890 >    pe = 0.0E0_DP
891 >    tauTemp = 0.0_dp
892 >    fTemp = 0.0_dp
893 >    tTemp = 0.0_dp
894   #else
895 <       tau = tauTemp
896 < #endif      
897 <    endif
895 >    qRow = 0.0_dp
896 >    qCol = 0.0_dp
897 >    
898 >    muRow = 0.0_dp
899 >    muCol = 0.0_dp
900 >    
901 >    u_lRow = 0.0_dp
902 >    u_lCol = 0.0_dp
903 >    
904 >    ARow = 0.0_dp
905 >    ACol = 0.0_dp
906 >    
907 >    fRow = 0.0_dp
908 >    fCol = 0.0_dp
909 >    
910 >  
911 >    tRow = 0.0_dp
912 >    tCol = 0.0_dp
913 >
914 >  
915  
916 +    eRow = 0.0_dp
917 +    eCol = 0.0_dp
918 +    eTemp = 0.0_dp
919 + #endif
920  
921 <  end subroutine do_force_loop
921 >  end subroutine zero_module_variables
922  
923 + #ifdef IS_MPI
924 + !! Function to properly build neighbor lists in MPI using newtons 3rd law.
925 + !! We don't want 2 processors doing the same i j pair twice.
926 + !! Also checks to see if i and j are the same particle.
927 +  function mpi_cycle_jLoop(i,j) result(do_cycle)
928 + !--------------- Arguments--------------------------
929 + ! Index i
930 +    integer,intent(in) :: i
931 + ! Index j
932 +    integer,intent(in) :: j
933 + ! Result do_cycle
934 +    logical :: do_cycle
935 + !--------------- Local variables--------------------
936 +    integer :: tag_i
937 +    integer :: tag_j
938 + !--------------- END DECLARATIONS------------------    
939 +    tag_i = tagRow(i)
940 +    tag_j = tagColumn(j)
941  
942 +    do_cycle = .false.
943  
944 +    if (tag_i == tag_j) then
945 +       do_cycle = .true.
946 +       return
947 +    end if
948 +
949 +    if (tag_i < tag_j) then
950 +       if (mod(tag_i + tag_j,2) == 0) do_cycle = .true.
951 +       return
952 +    else                
953 +       if (mod(tag_i + tag_j,2) == 1) do_cycle = .true.
954 +    endif
955 +  end function mpi_cycle_jLoop
956 + #endif
957 +
958   end module do_Forces

Diff Legend

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