ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/md_code/lj_FF.F90
(Generate patch)

Comparing trunk/mdtools/md_code/lj_FF.F90 (file contents):
Revision 280 by chuckv, Tue Feb 4 20:16:08 2003 UTC vs.
Revision 281 by chuckv, Mon Feb 24 21:25:15 2003 UTC

# Line 2 | Line 2
2   !! Corresponds to the force field defined in lj_FF.cpp
3   !! @author Charles F. Vardeman II
4   !! @author Matthew Meineke
5 < !! @version $Id: lj_FF.F90,v 1.18 2003-02-04 20:16:08 chuckv Exp $, $Date: 2003-02-04 20:16:08 $, $Name: not supported by cvs2svn $, $Revision: 1.18 $
5 > !! @version $Id: lj_FF.F90,v 1.19 2003-02-24 21:25:15 chuckv Exp $, $Date: 2003-02-24 21:25:15 $, $Name: not supported by cvs2svn $, $Revision: 1.19 $
6  
7  
8  
# Line 348 | Line 348 | contains
348  
349   !! FORCE routine Calculates Lennard Jones forces.
350   !------------------------------------------------------------->
351 <  subroutine do_lj_ff(q,f,potE,do_pot)
351 >  subroutine do_lj_ff(q,f,potE,tau,do_pot)
352   !! Position array provided by C, dimensioned by getNlocal
353      real ( kind = dp ), dimension(3,getNlocal()) :: q
354   !! Force array provided by C, dimensioned by getNlocal
355      real ( kind = dp ), dimension(3,getNlocal()) :: f
356 + !! Stress Tensor
357 +    real( kind = dp), dimension(9) :: tau
358 +    real( kind = dp), dimension(9) :: tauTemp
359      real ( kind = dp ) :: potE
360      logical ( kind = 2) :: do_pot
361      
362      type(lj_atype), pointer :: ljAtype_i
363      type(lj_atype), pointer :: ljAtype_j
364  
365 < #ifdef IS_MPI
366 <  real( kind = DP ), dimension(3,getNcol(plan_col)) :: efr
367 <  real( kind = DP ) :: pot_local
365 < #else
366 <  real( kind = DP ), dimension(3,getNlocal()) :: efr
367 < #endif
365 >
366 >
367 >
368    
369 < !! Local arrays needed for MPI
369 >
370   #ifdef IS_MPI
371 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow
372 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol
371 >  real( kind = DP ) :: pot_local
372  
373 <  real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow
374 <  real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol
375 <  real(kind = dp), dimension(3,getNlocal()) :: fMPITemp
373 > !! Local arrays needed for MPI
374 >  real(kind = dp), dimension(3,getNrow(plan_row)) :: qRow = 0.0_dp
375 >  real(kind = dp), dimension(3,getNcol(plan_col)) :: qCol = 0.0_dp
376 >
377 >  real(kind = dp), dimension(3,getNrow(plan_row)) :: fRow = 0.0_dp
378 >  real(kind = dp), dimension(3,getNcol(plan_col)) :: fCol = 0.0_dp
379 >  real(kind = dp), dimension(3,getNlocal()) :: fMPITemp = 0.0_dp
380  
381 <  real(kind = dp), dimension(getNrow(plan_row)) :: eRow
382 <  real(kind = dp), dimension(getNcol(plan_col)) :: eCol
381 >  real(kind = dp), dimension(getNrow(plan_row)) :: eRow = 0.0_dp
382 >  real(kind = dp), dimension(getNcol(plan_col)) :: eCol = 0.0_dp
383  
384 <  real(kind = dp), dimension(getNlocal()) :: eTemp
382 < #endif
384 >  real(kind = dp), dimension(getNlocal()) :: eTemp = 0.0_dp
385  
386 + #endif
387  
388  
389 +
390    real( kind = DP )   :: pe
391    logical             :: update_nlist
392  
# Line 392 | Line 396 | contains
396    integer :: j_start
397    integer :: tag_i,tag_j
398    real( kind = DP ) ::  r, pot, ftmp, dudr, d2, drdx1, kt1, kt2, kt3, ktmp
399 +  real( kind = dp ) :: fx,fy,fz
400 +  real( kind = DP ) ::  drdx, drdy, drdz
401    real( kind = DP ) ::  rxi, ryi, rzi, rxij, ryij, rzij, rijsq
402    real( kind = DP ) ::  rlistsq, rcutsq,rlist,rcut
403  
# Line 404 | Line 410 | contains
410    integer :: nrow
411    integer :: ncol
412    integer :: natoms
413 + !! should we calculate the stress tensor
414 +  logical  :: do_stress = .false.
415  
416  
417  
410 –
418   ! Make sure we are properly initialized.
419    if (.not. isljFFInit) then
420       write(default_error,*) "ERROR: lj_FF has not been properly initialized"
# Line 424 | Line 431 | contains
431    natoms = getNlocal()
432    call getRcut(rcut,rcut2=rcutsq)
433    call getRlist(rlist,rlistsq)
434 + !! Find ensemble
435 +  if (isEnsemble("NPT")) do_stress = .true.
436  
437   #ifndef IS_MPI
438    nrow = natoms - 1
# Line 462 | Line 471 | contains
471      eCol = 0.0E0_DP
472      eTemp = 0.0E0_DP
473   #endif
465 –    efr = 0.0E0_DP
474  
475   ! communicate MPI positions
476   #ifdef IS_MPI    
# Line 551 | Line 559 | contains
559                      pe = pe + pot
560   #endif                
561              
562 <                 efr(1,j) = -rxij
563 <                 efr(2,j) = -ryij
564 <                 efr(3,j) = -rzij
565 <
566 <                 do dim = 1, 3  
567 <
568 <            
569 <                    drdx1 = efr(dim,j) / r
562 <                    ftmp = dudr * drdx1
563 <
564 <
562 >                drdx = -rxij / r
563 >                drdy = -ryij / r
564 >                drdz = -rzij / r
565 >                
566 >                fx = dudr * drdx
567 >                fy = dudr * drdy
568 >                fz = dudr * drdz
569 >                
570   #ifdef IS_MPI
571 <                    fCol(dim,j) = fCol(dim,j) - ftmp
572 <                    fRow(dim,i) = fRow(dim,i) + ftmp
573 < #else                    
574 <            
575 <                    f(dim,j) = f(dim,j) - ftmp
576 <                    f(dim,i) = f(dim,i) + ftmp
577 <
578 < #endif                    
579 <                 enddo
580 <              endif
581 <           endif
582 <        enddo inner
571 >                fCol(1,j) = fCol(1,j) - fx
572 >                fCol(2,j) = fCol(2,j) - fy
573 >                fCol(3,j) = fCol(3,j) - fz
574 >                
575 >                fRow(1,j) = fRow(1,j) + fx
576 >                fRow(2,j) = fRow(2,j) + fy
577 >                fRow(3,j) = fRow(3,j) + fz
578 > #else
579 >                f(1,j) = f(1,j) - fx
580 >                f(2,j) = f(2,j) - fy
581 >                f(3,j) = f(3,j) - fz
582 >                f(1,i) = f(1,i) + fx
583 >                f(2,i) = f(2,i) + fy
584 >                f(3,i) = f(3,i) + fz
585 > #endif
586 >                
587 >                if (do_stress) then
588 >                   tauTemp(1) = tauTemp(1) + fx * rxij
589 >                   tauTemp(2) = tauTemp(2) + fx * ryij
590 >                   tauTemp(3) = tauTemp(3) + fx * rzij
591 >                   tauTemp(4) = tauTemp(4) + fy * rxij
592 >                   tauTemp(5) = tauTemp(5) + fy * ryij
593 >                   tauTemp(6) = tauTemp(6) + fy * rzij
594 >                   tauTemp(7) = tauTemp(7) + fz * rxij
595 >                   tauTemp(8) = tauTemp(8) + fz * ryij
596 >                   tauTemp(9) = tauTemp(9) + fz * rzij
597 >                endif
598 >             endif
599 >          enddo inner
600       enddo
601  
602   #ifdef IS_MPI
# Line 628 | Line 650 | contains
650   #else
651                  pe = pe + pot
652   #endif                
653 <
654 <                
655 <                 efr(1,j) = -rxij
656 <                 efr(2,j) = -ryij
657 <                 efr(3,j) = -rzij
658 <
659 <                 do dim = 1, 3                        
660 <                    
661 <                    drdx1 = efr(dim,j) / r
640 <                    ftmp = dudr * drdx1
653 >  
654 >                drdx = -rxij / r
655 >                drdy = -ryij / r
656 >                drdz = -rzij / r
657 >                
658 >                fx = dudr * drdx
659 >                fy = dudr * drdy
660 >                fz = dudr * drdz
661 >                
662   #ifdef IS_MPI
663 <                    fCol(dim,j) = fCol(dim,j) - ftmp
664 <                    fRow(dim,i) = fRow(dim,i) + ftmp
665 < #else                    
666 <                    f(dim,j) = f(dim,j) - ftmp
667 <                    f(dim,i) = f(dim,i) + ftmp
668 < #endif                    
669 <                 enddo
670 <              endif
671 <           enddo
672 <        endif
673 <     enddo
674 <  endif
663 >                fCol(1,j) = fCol(1,j) - fx
664 >                fCol(2,j) = fCol(2,j) - fy
665 >                fCol(3,j) = fCol(3,j) - fz
666 >                
667 >                fRow(1,j) = fRow(1,j) + fx
668 >                fRow(2,j) = fRow(2,j) + fy
669 >                fRow(3,j) = fRow(3,j) + fz
670 > #else
671 >                f(1,j) = f(1,j) - fx
672 >                f(2,j) = f(2,j) - fy
673 >                f(3,j) = f(3,j) - fz
674 >                f(1,i) = f(1,i) + fx
675 >                f(2,i) = f(2,i) + fy
676 >                f(3,i) = f(3,i) + fz
677 > #endif
678 >                
679 >                if (do_stress) then
680 >                   tauTemp(1) = tauTemp(1) + fx * rxij
681 >                   tauTemp(2) = tauTemp(2) + fx * ryij
682 >                   tauTemp(3) = tauTemp(3) + fx * rzij
683 >                   tauTemp(4) = tauTemp(4) + fy * rxij
684 >                   tauTemp(5) = tauTemp(5) + fy * ryij
685 >                   tauTemp(6) = tauTemp(6) + fy * rzij
686 >                   tauTemp(7) = tauTemp(7) + fz * rxij
687 >                   tauTemp(8) = tauTemp(8) + fz * ryij
688 >                   tauTemp(9) = tauTemp(9) + fz * rzij
689 >                endif
690 >                
691 >                
692 >             endif
693 >          enddo
694 >       endif
695 >    enddo
696 > endif
697 >
698  
699  
656 –
700   #ifdef IS_MPI
701      !!distribute forces
702  
# Line 668 | Line 711 | contains
711  
712      
713      if (do_pot) then
714 <    
672 < ! scatter/gather pot_row into the members of my column
714 >       ! scatter/gather pot_row into the members of my column
715         call scatter(eRow,eTemp,plan_row)
716        
717         ! scatter/gather pot_local into all other procs
# Line 684 | Line 726 | contains
726               pe_local = pe_local + eTemp(i)
727            enddo
728         endif
687 –    
729         pe = pe_local
730      endif
731 +
732   #endif
733  
734      potE = pe
735  
736 +
737 +    if (do_stress) then
738 + #ifdef IS_MPI
739 +       mpi_allreduce = (tau,tauTemp,9,mpi_double_precision,mpi_sum, &
740 +            mpi_comm_world,mpi_err)
741 + #else
742 +       tau = tauTemp
743 + #endif      
744 +    endif
745 +
746 +
747    end subroutine do_lj_ff
748  
749   !! Calculates the potential between two lj particles based on two lj_atype pointers, optionally returns second
# Line 786 | Line 839 | contains
839      real(kind = dp)  :: param1
840      real(kind = dp)  :: param2
841      real(kind = dp ) :: myMixParam
842 +    character(len = getStringLen()) :: thisMixingRule
843      integer, optional :: status
844  
845 <
845 > !! get the mixing rules from the simulation
846 >    thisMixingRule = returnMixingRules()
847      myMixParam = 0.0_dp
848  
849      if (present(status)) status = 0
850 <
851 <    select case (thisParam)
852 <
853 <    case ("sigma")
854 <       myMixParam = 0.5_dp * (param1 + param2)
855 <    case ("epsilon")
856 <       myMixParam = sqrt(param1 * param2)
850 >    select case (thisMixingRule)
851 >       case ("standard")
852 >          select case (thisParam)
853 >          case ("sigma")
854 >             myMixParam = 0.5_dp * (param1 + param2)
855 >          case ("epsilon")
856 >             myMixParam = sqrt(param1 * param2)
857 >          case default
858 >             status = -1
859 >          end select
860 >    case("LJglass")
861      case default
862         status = -1
863      end select
805 –
864    end function calcLJMix
865  
866  

Diff Legend

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