FIST nonbonded neighbor lists (modified version)

svn-origin-rev: 991
This commit is contained in:
Matthias Krack 2002-11-21 17:13:38 +00:00
parent 2bd9e96ec4
commit dbf2f16008
10 changed files with 1321 additions and 73 deletions

View file

@ -85,6 +85,7 @@ OBJECTS_GENERIC =\
fist_force_numer.o\
fist_intra_force.o\
fist_main.o\
fist_neighbor_lists.o\
fist_nonbond_force.o\
force_control.o\
force_fields.o\

View file

@ -218,7 +218,9 @@ SUBROUTINE force_control ( rep_env, ewald_param, potparm, thermo, &
!
f_nonbond = zero
CALL force_nonbond ( ewald_param, part, pnode, box, potparm, &
pot_nonbond, f_nonbond, pv_nonbond )
pot_nonbond, f_nonbond, pv_nonbond, &
rep_env%ll_data(1)%nonbonded, &
rep_env%ll_data(1)%r_last_update )
!
! get g-space non-bonded forces:
!

View file

@ -312,6 +312,7 @@ SUBROUTINE fist ( globenv )
rep_env ( ibead ) % box % unit_of_length_name = "ANGSTROM"
rep_env ( ibead ) % box % unit_of_length = 1.0_dbl
rep_env ( ibead ) % box % scaled_coordinates = .FALSE.
rep_env ( ibead ) % box % subcells = simpar % subcells
! dbg TEMPORARY FIX
! If run is a debug. Make box_ref = box to work under

1148
src/fist_neighbor_lists.F Normal file

File diff suppressed because it is too large Load diff

View file

@ -30,17 +30,23 @@ MODULE fist_nonbond_force
linklist_atoms, linklist_exclusion
USE pair_potential, ONLY : potential_s, potentialparm_type
USE particle_types, ONLY : particle_type
USE qs_neighbor_list_types, ONLY: first_node,&
USE qs_neighbor_list_types, ONLY: first_list,&
first_node,&
get_neighbor_list,&
get_neighbor_list_set,&
get_neighbor_node,&
neighbor_list_set_type,&
neighbor_list_type,&
neighbor_node_type,&
next
USE simulation_cell, ONLY : cell_type, pbc, get_cell_param, &
real_to_scaled, scaled_to_real
real_to_scaled, scaled_to_real,pbc
USE termination, ONLY : stop_memory
USE timings, ONLY : timeset, timestop
USE util, ONLY : include_list
USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
IMPLICIT NONE
PRIVATE
@ -53,7 +59,8 @@ CONTAINS
!******************************************************************************
SUBROUTINE force_nonbond ( ewald_param, part, pnode, box, potparm, &
pot_nonbond, f_nonbond, ptens_nonbond )
pot_nonbond, f_nonbond, ptens_nonbond, &
nonbonded, r_last_update )
! Calculates the force and the potential of the minimum image, and
! the pressure tensor
@ -70,19 +77,28 @@ SUBROUTINE force_nonbond ( ewald_param, part, pnode, box, potparm, &
REAL ( dbl ), INTENT ( OUT ), DIMENSION ( :, : ) :: f_nonbond
REAL ( dbl ), INTENT ( OUT ), DIMENSION ( :, : ) :: ptens_nonbond
TYPE(neighbor_list_set_p_type), DIMENSION(:), INTENT(IN) :: nonbonded
REAL(dbl), DIMENSION(:,:), INTENT(IN) :: r_last_update
! Locals
CHARACTER(LEN=*), PARAMETER :: routine_name = "force_nonbond"
CHARACTER(LEN=*), PARAMETER :: module_name = "fist_nonbond_force"
CHARACTER(LEN=*), PARAMETER :: routine =&
"SUBROUTINE "//routine_name//" (MODULE "//module_name//")"
INTEGER :: ikind, jkind, nkinds, inode, natoms, nnodes
INTEGER :: ikind, jkind, nkinds, inode, natoms, nnodes, ab
INTEGER :: handle, istat, nneighbor, ineighbor, atom_a, atom_b, cell_b ( 3 )
REAL ( dbl ), DIMENSION (3) :: sab_pbc, rab, rb, sab, ra
REAL ( dbl ) :: energy, fscalar, rab2, flops
REAL ( dbl ) :: energy, fscalar, rab2, flops, rab2_max
REAL ( dbl ), ALLOCATABLE, DIMENSION ( :, : ) :: rtest
TYPE ( neighbor_list_type ), POINTER :: neighbor_list
TYPE ( neighbor_node_type ), POINTER :: neighbor_node
TYPE ( cell_type ), POINTER :: cell
!MK
TYPE(neighbor_list_set_type), POINTER :: neighbor_list_set
INTEGER :: ilist,inode,nlist,nnode
REAL(dbl), DIMENSION(3) :: dra,drb,rab_last_update
!------------------------------------------------------------------------------
@ -97,55 +113,111 @@ SUBROUTINE force_nonbond ( ewald_param, part, pnode, box, potparm, &
! local copy of cutoffs
nkinds = SIZE ( potparm, 1 )
ALLOCATE ( rtest ( nkinds, nkinds ), STAT = istat )
IF ( istat /= 0 ) CALL stop_memory ( routine, "rtest", nkinds ** 2 )
DO ikind = 1, nkinds
DO jkind = 1, nkinds
rtest ( ikind, jkind ) = potparm ( ikind, jkind ) % rcutsq
END DO
END DO
!MK ALLOCATE ( rtest ( nkinds, nkinds ), STAT = istat )
!MK IF ( istat /= 0 ) CALL stop_memory ( routine, "rtest", nkinds ** 2 )
!MK
!MK DO ikind = 1, nkinds
!MK DO jkind = 1, nkinds
!MK rtest ( ikind, jkind ) = potparm ( ikind, jkind ) % rcutsq
!MK END DO
!MK END DO
!
! starting the force loop
!
nnodes = SIZE ( pnode )
DO inode = 1, nnodes
atom_a = pnode ( inode ) % p % iatom
neighbor_list => pnode ( inode ) % nl
nneighbor = pnode ( inode ) % nneighbor
ra ( : ) = part ( atom_a ) % r ( : )
CALL get_atomic_kind ( part ( atom_a ) % atomic_kind, &
kind_number = ikind )
! now do neighbors
neighbor_node => first_node ( neighbor_list )
DO ineighbor = 1, nneighbor
CALL get_neighbor_node ( neighbor_node = neighbor_node, &
neighbor = atom_b, &
cell = cell_b )
rb ( : ) = part ( atom_b ) % r ( : )
CALL get_atomic_kind ( part ( atom_b ) % atomic_kind, &
kind_number = jkind )
rab ( : ) = rb ( : ) - ra ( : )
sab ( : ) = real_to_scaled(rab(:),cell)
sab_pbc ( : ) = sab ( : ) + REAL ( cell_b ( : ), dbl )
rab ( : ) = scaled_to_real(sab_pbc(:),cell)
rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
!MK
DO ikind=1,nkinds
DO jkind=ikind,nkinds
IF ( rab2 <= rtest ( ikind, jkind ) ) THEN
CALL potential_s ( rab2, ikind, jkind, energy, fscalar )
!
! summing up the potential energy,the force and pressure tensor
!
CALL sum_ener_forces ( atom_a, atom_b, pot_nonbond, energy, &
fscalar, f_nonbond, rab, ptens_nonbond )
flops = flops + 64.0_dbl
END IF
ab = ikind + jkind*(jkind - 1)/2
neighbor_node => next ( neighbor_node )
neighbor_list_set => nonbonded(ab)%neighbor_list_set
END DO
END DO
IF (.NOT.ASSOCIATED(neighbor_list_set)) CYCLE
rab2_max = potparm(ikind,jkind)%rcutsq
CALL get_neighbor_list_set(neighbor_list_set=neighbor_list_set,&
nlist=nlist)
neighbor_list => first_list(neighbor_list_set)
DO ilist=1,nlist
CALL get_neighbor_list(neighbor_list=neighbor_list,&
atom=atom_a,&
nnode=nnode)
dra(:) = part(atom_a)%r(:) - r_last_update(:,atom_a)
neighbor_node => first_node(neighbor_list)
DO inode=1,nnode
CALL get_neighbor_node(neighbor_node=neighbor_node,&
neighbor=atom_b,&
cell=cell_b,&
r=rab_last_update)
drb(:) = part(atom_b)%r(:) - r_last_update(:,atom_b)
rab(:) = rab_last_update(:) - dra(:) + drb(:)
rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
IF (rab2 <= rab2_max) THEN
CALL potential_s(rab2,ikind,jkind,energy,fscalar)
CALL sum_ener_forces(atom_a,atom_b,pot_nonbond,energy,&
fscalar,f_nonbond,rab,ptens_nonbond)
END IF
neighbor_node => next(neighbor_node)
END DO
neighbor_list => next(neighbor_list)
END DO
END DO
END DO
!MK nnodes = SIZE ( pnode )
!MK DO inode = 1, nnodes
!MK atom_a = pnode ( inode ) % p % iatom
!MK neighbor_list => pnode ( inode ) % nl
!MK nneighbor = pnode ( inode ) % nneighbor
!MK ra ( : ) = part ( atom_a ) % r ( : )
!MK CALL get_atomic_kind ( part ( atom_a ) % atomic_kind, &
!MK kind_number = ikind )
!MK! now do neighbors
!MK neighbor_node => first_node ( neighbor_list )
!MK DO ineighbor = 1, nneighbor
!MK CALL get_neighbor_node ( neighbor_node = neighbor_node, &
!MK neighbor = atom_b, &
!MK cell = cell_b )
!MK rb ( : ) = part ( atom_b ) % r ( : )
!MK CALL get_atomic_kind ( part ( atom_b ) % atomic_kind, &
!MK kind_number = jkind )
!MK rab ( : ) = rb ( : ) - ra ( : )
!MK sab ( : ) = real_to_scaled(rab(:),cell)
!MK sab_pbc ( : ) = sab ( : ) + REAL ( cell_b ( : ), dbl )
!MK rab ( : ) = scaled_to_real(sab_pbc(:),cell)
!MK rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
!MK
!MK IF ( rab2 <= rtest ( ikind, jkind ) ) THEN
!MK CALL potential_s ( rab2, ikind, jkind, energy, fscalar )
!MK!
!MK! summing up the potential energy,the force and pressure tensor
!MK!
!MK CALL sum_ener_forces ( atom_a, atom_b, pot_nonbond, energy, &
!MK fscalar, f_nonbond, rab, ptens_nonbond )
!MK flops = flops + 64.0_dbl
!MK END IF
!MK
!MK neighbor_node => next ( neighbor_node )
!MK
!MK END DO
!MK END DO
! computing long range corrections to the potential
!
@ -155,10 +227,10 @@ SUBROUTINE force_nonbond ( ewald_param, part, pnode, box, potparm, &
! pot_nonbond=pot_nonbond+lrc*(1./box % deth)
!
DEALLOCATE ( rtest, STAT = istat )
IF ( istat /= 0 ) CALL stop_memory ( routine, "rtest" )
!
flops = flops * 1.E-6_dbl
!MK DEALLOCATE ( rtest, STAT = istat )
!MK IF ( istat /= 0 ) CALL stop_memory ( routine, "rtest" )
!MK
!MK flops = flops * 1.E-6_dbl
CALL timestop ( flops, handle )
END SUBROUTINE force_nonbond

View file

@ -33,6 +33,10 @@ MODULE linklist_control
USE termination, ONLY : stop_memory, stop_program
USE timings, ONLY : timeset, timestop
USE fist_neighbor_lists, ONLY: build_fist_neighbor_lists,&
rebuild_fist_neighbor_lists
USE pair_potential, ONLY: potentialparm_type
IMPLICIT NONE
PRIVATE
@ -81,7 +85,8 @@ SUBROUTINE list_control ( rep_env )
internal_data => rep_env % ll_data ( 1 )
atomic_kind_set => rep_env % atomic_kind_set
first_time = .NOT.ASSOCIATED ( internal_data % sab_ppnl )
!MK first_time = .NOT.ASSOCIATED ( internal_data % sab_ppnl )
first_time = .NOT.ASSOCIATED ( internal_data % nonbonded )
nnodes = SIZE ( pnode )
IF ( .NOT. ASSOCIATED ( internal_data % r_last_update ) ) THEN
@ -115,13 +120,20 @@ SUBROUTINE list_control ( rep_env )
list_update_flag = .FALSE.
END IF
IF ( list_update_flag ) THEN
IF ( first_time ) THEN
CALL build_verlet_lists (atomic_kind_set, part, box, internal_data, pnode )
CALL build_fist_neighbor_lists(atomic_kind_set,part,pnode,box,&
internal_data%rlist_cut,&
internal_data%nonbonded,&
internal_data%globenv)
!MK CALL build_verlet_lists (atomic_kind_set, part, box, internal_data, pnode )
ELSE
CALL rebuild_verlet_lists ( atomic_kind_set, part, box, internal_data, pnode )
CALL rebuild_fist_neighbor_lists(atomic_kind_set,part,pnode,box,&
internal_data%rlist_cut,&
internal_data%nonbonded,&
internal_data%globenv)
!MK CALL rebuild_verlet_lists ( atomic_kind_set, part, box, internal_data, pnode )
ENDIF
iw = internal_data % globenv % scr

View file

@ -41,16 +41,16 @@ MODULE linklist_types
END TYPE subcell_data_type
TYPE linklist_internal_data_type
TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: nonbonded
TYPE ( neighbor_list_set_p_type ), DIMENSION ( : ), POINTER :: sab_ppnl
TYPE ( subcell_data_type ), DIMENSION ( :, : ), POINTER :: ll_cell
TYPE ( global_environment_type ) :: globenv
INTEGER :: natom_types
INTEGER :: counter, last_update, num_update
INTEGER :: print_level
INTEGER :: subcells
REAL ( dbl ), DIMENSION ( :, : ), POINTER :: rlist_cut
REAL ( dbl ), DIMENSION ( :, : ), POINTER :: rlist_cutsq
REAL ( dbl ) :: verlet_skin
REAL ( dbl ) :: subcells, verlet_skin
REAL ( dbl ) :: lup, aup
REAL ( dbl ), DIMENSION ( :, : ), POINTER :: r_last_update
END TYPE linklist_internal_data_type
@ -68,9 +68,8 @@ SUBROUTINE set_linklist_internal_data ( internal_data, globenv, vskin, &
! Arguments
TYPE ( linklist_internal_data_type ), INTENT ( INOUT ) :: internal_data
TYPE ( global_environment_type ), INTENT ( IN ) :: globenv
REAL ( dbl ), INTENT ( IN ), OPTIONAL :: vskin
REAL ( dbl ), INTENT ( IN ), OPTIONAL :: subcells, vskin
REAL ( dbl ), INTENT ( IN ), OPTIONAL :: rcut ( :, : )
INTEGER, INTENT ( IN ), OPTIONAL :: subcells
INTEGER, INTENT ( IN ), OPTIONAL :: natype
INTEGER, INTENT ( IN ), OPTIONAL :: count
INTEGER, INTENT ( IN ), OPTIONAL :: printlevel
@ -155,6 +154,7 @@ SUBROUTINE initialize_linklist_data ( internal_data )
NULLIFY ( internal_data % rlist_cut )
NULLIFY ( internal_data % rlist_cutsq )
NULLIFY ( internal_data % sab_ppnl )
NULLIFY ( internal_data % nonbonded )
NULLIFY ( internal_data % ll_cell )
END SUBROUTINE initialize_linklist_data

View file

@ -109,9 +109,9 @@ CONTAINS
LOGICAL :: ionode
INTEGER :: group, mype, npe, ikind, inode, &
nkind, nnodes, istat, handle, &
subcells, atom_a, natoms, periodic ( 3 )
atom_a, natoms, periodic ( 3 )
INTEGER, DIMENSION ( : ), ALLOCATABLE :: natom
REAL ( dbl ) :: ra ( 3 ), lat_vec ( 3 ), sa ( 3 ), unit_of_length
REAL ( dbl ) :: ra ( 3 ), lat_vec ( 3 ), sa ( 3 ), unit_of_length, subcells
TYPE ( atomic_kind_type ), POINTER :: atomic_kind
TYPE ( atoms_type ), DIMENSION ( : ), ALLOCATABLE :: atoms
TYPE ( global_environment_type ) :: globenv
@ -232,7 +232,7 @@ CONTAINS
! History: - Creation (20.03.2002,MK)
! ***************************************************************************
INTEGER, INTENT ( IN ) :: subcells
REAL ( dbl ), INTENT ( IN ) :: subcells
REAL ( dbl ), INTENT ( IN ) :: rlist_cutsq ( :, : )
REAL ( dbl ), INTENT ( IN ) :: rlist_cut ( :, : )
INTEGER, INTENT ( IN ) :: natom ( : )
@ -618,9 +618,9 @@ CONTAINS
LOGICAL :: ionode
INTEGER :: group, mype, npe, ikind, inode, &
nkind, nnodes, istat, handle, output_unit, &
subcells, atom_a, natoms, periodic ( 3 )
atom_a, natoms, periodic ( 3 )
INTEGER, DIMENSION ( : ), ALLOCATABLE :: natom
REAL ( dbl ) :: ra ( 3 ), sa ( 3 ), lat_vec ( 3 )
REAL ( dbl ) :: ra ( 3 ), sa ( 3 ), lat_vec ( 3 ), subcells
TYPE ( atomic_kind_type ), POINTER :: atomic_kind
TYPE ( atoms_type ), DIMENSION ( : ), ALLOCATABLE :: atoms
TYPE ( global_environment_type ) :: globenv

View file

@ -62,7 +62,7 @@ MODULE md
INTEGER :: nc
INTEGER :: nyosh
INTEGER :: nhclen
INTEGER :: subcells
REAL ( dbl ) :: subcells
REAL ( dbl ) :: tau_nhc
REAL ( dbl ) :: tau_cell
REAL ( dbl ), POINTER, DIMENSION ( : ) :: dt_yosh
@ -173,7 +173,7 @@ SUBROUTINE read_md_section ( simpar, mdpar, mdio )
simpar % shake_tol = 1.0E-6_dbl
simpar % nhclen = 1
simpar % nc = 1
simpar % subcells = 0
simpar % subcells = 2.0_dbl
simpar % nyosh = 1
simpar % tau_nhc = 1000.0_dbl
simpar % tau_cell = 1000.0_dbl

View file

@ -219,7 +219,7 @@ CONTAINS
! *****************************************************************************
SUBROUTINE add_neighbor_node(neighbor_list,neighbor,r,cell)
SUBROUTINE add_neighbor_node(neighbor_list,neighbor,r,cell,exclusion_list)
! Purpose: Add a new neighbor list node to neighbor list.
@ -227,10 +227,11 @@ CONTAINS
! ***************************************************************************
TYPE(neighbor_list_type), POINTER :: neighbor_list
INTEGER, INTENT(IN) :: neighbor
REAL(wp), DIMENSION(3), INTENT(IN) :: r
INTEGER, DIMENSION(3), INTENT(IN) :: cell
TYPE(neighbor_list_type), POINTER :: neighbor_list
INTEGER, INTENT(IN) :: neighbor
REAL(wp), DIMENSION(3), INTENT(IN) :: r
INTEGER, DIMENSION(3), INTENT(IN) :: cell
INTEGER, DIMENSION(:), OPTIONAL, INTENT(IN) :: exclusion_list
! *** Local parameters ***
@ -240,10 +241,19 @@ CONTAINS
! *** Local variables ***
TYPE(neighbor_node_type), POINTER :: new_neighbor_node
INTEGER :: istat
INTEGER :: iatom,istat
! ---------------------------------------------------------------------------
! *** Check for exclusions ***
IF (PRESENT(exclusion_list)) THEN
DO iatom=1,SIZE(exclusion_list)
IF (exclusion_list(iatom) == 0) EXIT
IF (exclusion_list(iatom) == neighbor) RETURN
END DO
END IF
IF (ASSOCIATED(neighbor_list%last_neighbor_node)) THEN
new_neighbor_node => neighbor_list%last_neighbor_node%next_neighbor_node
@ -374,6 +384,7 @@ CONTAINS
IF (ASSOCIATED(neighbor_list_set%last_neighbor_list)) THEN
neighbor_list => neighbor_list_set%last_neighbor_list%next_neighbor_list
NULLIFY (neighbor_list_set%last_neighbor_list%next_neighbor_list)
ELSE
neighbor_list => neighbor_list_set%first_neighbor_list
NULLIFY (neighbor_list_set%first_neighbor_list)
@ -395,6 +406,7 @@ CONTAINS
IF (ASSOCIATED(neighbor_list%last_neighbor_node)) THEN
neighbor_node => neighbor_list%last_neighbor_node%next_neighbor_node
NULLIFY (neighbor_list%last_neighbor_node%next_neighbor_node)
ELSE
neighbor_node => neighbor_list%first_neighbor_node
NULLIFY (neighbor_list%first_neighbor_node)