1 |
|
!! Calculates Long Range forces. |
2 |
|
!! @author Charles F. Vardeman II |
3 |
|
!! @author Matthew Meineke |
4 |
< |
!! @version $Id: do_Forces.F90,v 1.3 2003-03-06 19:57:03 chuckv Exp $, $Date: 2003-03-06 19:57:03 $, $Name: not supported by cvs2svn $, $Revision: 1.3 $ |
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 |
|
|
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(3,getNrow(plan_row)) :: ARow = 0.0_dp |
66 |
< |
real(kind = dp), dimension(3,getNcol(plan_col)) :: ACol = 0.0_dp |
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 |
– |
|
69 |
– |
|
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 |
|
|
73 |
– |
|
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 |
76 |
– |
|
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 |
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,getNlocal()) :: fTemp = 0.0_dp |
91 |
< |
real(kind = dp), dimension(3,getNlocal()) :: tTemp = 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. |
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. |
399 |
|
if (isPBC()) wrap = .true. |
400 |
|
|
401 |
|
|
396 |
– |
#ifndef IS_MPI |
397 |
– |
nrow = natoms - 1 |
398 |
– |
ncol = natoms |
399 |
– |
#else |
400 |
– |
nrow = getNrow(plan_row) |
401 |
– |
ncol = getNcol(plan_col) |
402 |
– |
nlocal = natoms |
403 |
– |
j_start = 1 |
404 |
– |
#endif |
402 |
|
|
403 |
|
|
404 |
|
!! See if we need to update neighbor lists |
428 |
|
#endif |
429 |
|
|
430 |
|
|
434 |
– |
if (update_nlist) then |
435 |
– |
|
436 |
– |
! save current configuration, contruct neighbor list, |
437 |
– |
! and calculate forces |
438 |
– |
call save_neighborList(q) |
439 |
– |
|
440 |
– |
neighborListSize = getNeighborListSize() |
441 |
– |
nlist = 0 |
442 |
– |
|
443 |
– |
|
444 |
– |
|
445 |
– |
do i = 1, nrow |
446 |
– |
point(i) = nlist + 1 |
431 |
|
#ifdef IS_MPI |
432 |
< |
Atype_i => identPtrListRow(i)%this |
433 |
< |
tag_i = tagRow(i) |
450 |
< |
#else |
451 |
< |
Atype_i => identPtrList(i)%this |
452 |
< |
j_start = i + 1 |
453 |
< |
#endif |
454 |
< |
|
455 |
< |
inner: do j = j_start, ncol |
456 |
< |
! Assign identity pointers and tags |
457 |
< |
#ifdef IS_MPI |
458 |
< |
Atype_j => identPtrListColumn(j)%this |
459 |
< |
|
460 |
< |
call get_interatomic_vector(i,j,qRow(:,i),qCol(:,j),& |
461 |
< |
rxij,ryij,rzij,rijsq,r) |
462 |
< |
!! For mpi, use newtons 3rd law when building neigbor list |
463 |
< |
!! Also check to see the particle i != j. |
464 |
< |
if (mpi_cycle_jLoop(i,j)) cycle inner: |
465 |
< |
|
466 |
< |
#else |
467 |
< |
Atype_j => identPtrList(j)%this |
468 |
< |
call get_interatomic_vector(i,j,q(:,i),q(:,j),& |
469 |
< |
rxij,ryij,rzij,rijsq,r) |
432 |
> |
|
433 |
> |
if (update_nlist) then |
434 |
|
|
435 |
< |
#endif |
436 |
< |
|
437 |
< |
if (rijsq < rlistsq) then |
438 |
< |
|
439 |
< |
nlist = nlist + 1 |
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 |
< |
if (nlist > neighborListSize) then |
483 |
< |
call expandList(listerror) |
484 |
< |
if (listerror /= 0) then |
480 |
< |
FFerror = -1 |
481 |
< |
write(DEFAULT_ERROR,*) "ERROR: nlist > list size and max allocations exceeded." |
482 |
< |
return |
483 |
< |
end if |
484 |
< |
endif |
482 |
> |
point(nrow + 1) = nlist + 1 |
483 |
> |
|
484 |
> |
else !! (update) |
485 |
|
|
486 |
< |
list(nlist) = j |
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 |
|
|
508 |
< |
if (rijsq < rcutsq) then |
509 |
< |
call do_pair(Atype_i,Atype_j,i,j,r,rxij,ryij,rzij) |
510 |
< |
endif |
511 |
< |
enddo inner |
512 |
< |
enddo |
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 |
< |
#ifdef IS_MPI |
523 |
< |
point(nrow + 1) = nlist + 1 |
524 |
< |
#else |
525 |
< |
point(natoms) = nlist + 1 |
526 |
< |
#endif |
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 |
< |
else !! (update) |
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 |
< |
! use the list to find the neighbors |
543 |
< |
do i = 1, nrow |
544 |
< |
JBEG = POINT(i) |
545 |
< |
JEND = POINT(i+1) - 1 |
546 |
< |
! check thiat molecule i has neighbors |
547 |
< |
if (jbeg .le. jend) then |
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 |
< |
ljAtype_i => identPtrListRow(i)%this |
578 |
< |
#else |
579 |
< |
ljAtype_i => identPtrList(i)%this |
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 |
< |
do jnab = jbeg, jend |
585 |
< |
j = list(jnab) |
584 |
> |
|
585 |
> |
! This is the post-pair loop: |
586 |
|
#ifdef IS_MPI |
587 |
< |
ljAtype_j = identPtrListColumn(j)%this |
588 |
< |
call get_interatomic_vector(i,j,qRow(:,i),qCol(:,j),& |
589 |
< |
rxij,ryij,rzij,rijsq,r) |
590 |
< |
|
587 |
> |
|
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 |
< |
ljAtype_j = identPtrList(j)%this |
597 |
< |
call get_interatomic_vector(i,j,q(:,i),q(:,j),& |
598 |
< |
rxij,ryij,rzij,rijsq,r) |
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 |
527 |
– |
call do_pair(i,j,r,rxij,ryij,rzij) |
528 |
– |
enddo |
529 |
– |
endif |
530 |
– |
enddo |
531 |
– |
endif |
532 |
– |
|
605 |
|
|
606 |
|
|
607 |
+ |
|
608 |
+ |
|
609 |
|
#ifdef IS_MPI |
610 |
|
!!distribute forces |
611 |
|
|
612 |
< |
call scatter(fRow,f,plan_row3d) |
613 |
< |
call scatter(fCol,fTemp,plan_col3d) |
612 |
> |
call scatter(fRow,fTemp1,plan_row3d) |
613 |
> |
call scatter(fCol,fTemp2,plan_col3d) |
614 |
|
|
615 |
+ |
|
616 |
|
do i = 1,nlocal |
617 |
< |
f(1:3,i) = f(1:3,i) + fTemp(1:3,i) |
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,t,plan_row3d) |
622 |
< |
call scatter(tCol,tTemp,plan_col3d) |
621 |
> |
call scatter(tRow,tTemp1,plan_row3d) |
622 |
> |
call scatter(tCol,tTemp2,plan_col3d) |
623 |
|
|
624 |
|
do i = 1,nlocal |
625 |
< |
t(1:3,i) = t(1:3,i) + tTemp(1:3,i) |
625 |
> |
tTemp(1:3,i) = tTemp1(1:3,i) + tTemp2(1:3,i) |
626 |
|
end do |
627 |
|
endif |
628 |
|
|
646 |
|
endif |
647 |
|
#else |
648 |
|
! Copy local array into return array for c |
649 |
< |
f = fTemp |
650 |
< |
t = tTemp |
649 |
> |
f = f+fTemp |
650 |
> |
t = t+tTemp |
651 |
|
#endif |
652 |
|
|
653 |
|
potE = pe |
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 |
< |
|
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 |
|
|
737 |
|
|
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 |
+ |
if (Atype_i%is_dp .and. Atype_j%is_dp) then |
743 |
|
|
744 |
+ |
#ifdef IS_MPI |
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 |
+ |
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 |
|
|
752 |
+ |
if (do_reaction_field) then |
753 |
+ |
#ifdef IS_MPI |
754 |
+ |
call accumulate_rf(i, j, r_ij, rflRow(:,i), rflCol(:j), & |
755 |
+ |
ulRow(:i), ulCol(:,j), rt, rrf) |
756 |
+ |
#else |
757 |
+ |
call accumulate_rf(i, j, r_ij, rfl(:,i), rfl(:j), & |
758 |
+ |
ul(:,i), ul(:,j), rt, rrf) |
759 |
+ |
#endif |
760 |
+ |
endif |
761 |
|
|
762 |
< |
call getLJForce(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 |