From dbf2f16008dc962cbfe94698812fd47d684a3c09 Mon Sep 17 00:00:00 2001 From: Matthias Krack Date: Thu, 21 Nov 2002 17:13:38 +0000 Subject: [PATCH] FIST nonbonded neighbor lists (modified version) svn-origin-rev: 991 --- src/OBJECTDEFS | 1 + src/fist_force.F | 4 +- src/fist_main.F | 1 + src/fist_neighbor_lists.F | 1148 ++++++++++++++++++++++++++++++++++ src/fist_nonbond_force.F | 174 ++++-- src/linklist_control.F | 20 +- src/linklist_types.F | 8 +- src/linklist_verlet_list.F | 10 +- src/md.F | 4 +- src/qs_neighbor_list_types.F | 24 +- 10 files changed, 1321 insertions(+), 73 deletions(-) create mode 100644 src/fist_neighbor_lists.F diff --git a/src/OBJECTDEFS b/src/OBJECTDEFS index 9c91108885..f094b161c3 100644 --- a/src/OBJECTDEFS +++ b/src/OBJECTDEFS @@ -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\ diff --git a/src/fist_force.F b/src/fist_force.F index e36d87a265..a326ec0778 100644 --- a/src/fist_force.F +++ b/src/fist_force.F @@ -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: ! diff --git a/src/fist_main.F b/src/fist_main.F index f618c8fb83..913a2bb539 100644 --- a/src/fist_main.F +++ b/src/fist_main.F @@ -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 diff --git a/src/fist_neighbor_lists.F b/src/fist_neighbor_lists.F new file mode 100644 index 0000000000..5b8ed31b94 --- /dev/null +++ b/src/fist_neighbor_lists.F @@ -0,0 +1,1148 @@ +!-----------------------------------------------------------------------------! +! CP2K: A general program to perform molecular dynamics simulations +! Copyright (C) 2001 - 2002 CP2K developers group +!-----------------------------------------------------------------------------! +!!****** cp2k/fist_neighbor_lists [1.0] * +!! +!! NAME +!! fist_neighbor_lists +!! +!! FUNCTION +!! Build all neighbor lists for FIST. +!! +!! AUTHOR +!! Matthias Krack (19.11.2002) +!! +!! MODIFICATION HISTORY +!! none +!! +!! SOURCE +!****************************************************************************** + +MODULE fist_neighbor_lists + + USE kinds, ONLY: int_size,& + wp => dp,& + wp_size => dp_size + + USE atomic_kind_types, ONLY: atomic_kind_type,& + get_atomic_kind,& + get_atomic_kind_set + USE checkpoint_handler, ONLY: write_checkpoint_information + USE global_types, ONLY: LOW,global_environment_type + USE molecule_types, ONLY: particle_node_type,& + linklist_exclusion + USE pair_potential, ONLY: potentialparm_type + USE particle_types, ONLY: particle_type + USE qs_neighbor_list_types, ONLY: add_neighbor_list,& + add_neighbor_node,& + allocate_neighbor_list_set,& + clean_neighbor_list_set,& + deallocate_neighbor_list_set,& + find_neighbor_list,& + first_list,& + first_node,& + get_neighbor_list,& + get_neighbor_node,& + init_neighbor_list,& + init_neighbor_list_set,& + neighbor_list_p_type,& + neighbor_list_set_p_type,& + neighbor_list_set_type,& + neighbor_list_type,& + neighbor_node_type,& + next,& + set_neighbor_node + USE simulation_cell, ONLY: cell_type,& + get_cell,& + pbc,& + real_to_scaled,& + scaled_to_real + USE termination, ONLY: stop_memory,& + stop_program + USE timings, ONLY: timeset,& + timestop + + IMPLICIT NONE + + PRIVATE + +! *** Global types of the module *** + + TYPE atoms_type + INTEGER, DIMENSION(:), POINTER :: list + END TYPE atoms_type + + TYPE excl_type + INTEGER, DIMENSION(:,:), POINTER :: list + END TYPE excl_type + + TYPE pbc_coord_type + REAL(wp), DIMENSION(:,:), POINTER :: r,s + END TYPE pbc_coord_type + +! *** Global parameters *** + + CHARACTER(LEN=*), PARAMETER :: module_name = "fist_neighbor_lists" + +! *** Global variables of the module *** + + TYPE(linklist_exclusion), POINTER :: excl_node + TYPE(neighbor_list_set_type), POINTER :: neighbor_list_set + TYPE(neighbor_list_type), POINTER :: neighbor_list + TYPE(neighbor_node_type), POINTER :: neighbor_node + CHARACTER(LEN=40) :: string + CHARACTER(LEN=8) :: unit_of_length_name + REAL(wp) :: a_max,a_min,& + b_max,b_min,& + c_max,c_min,& + rab_max,rab2,rab2_max,& + subcells,unit_of_length + INTEGER :: atom_a,atom_b,group,& + iatom,icell,iexcl,igrid,iijk,ikind,& + ineighbor,inode,ipe,istat,& + jatom,jcell,jgrid,jkind,& + kcell,kgrid,maxatom,maxexcl,mype,& + nkind,nneighbor,nnode,npe,& + output_unit + LOGICAL :: ionode,print_cell_parameters + + REAL(wp), DIMENSION(3) :: r,ra_pbc,rab,rb,s,s_max,s_min,sab,sab_max,sb,& + sb_max,sb_min,sb_pbc,sa_pbc + INTEGER, DIMENSION(3) :: cell_a,cell_b,ncell,ngrid,periodic + + TYPE(atoms_type), DIMENSION(:), ALLOCATABLE :: atoms + TYPE(excl_type), DIMENSION(:), ALLOCATABLE :: excl_a + TYPE(neighbor_list_p_type), DIMENSION(:), ALLOCATABLE :: kind_a + TYPE(pbc_coord_type), DIMENSION(:), ALLOCATABLE :: pbc_coord + REAL(wp), DIMENSION(:), ALLOCATABLE :: kind_radius + INTEGER, DIMENSION(:), ALLOCATABLE :: atom_of_kind,& + kind_of,& + natom,neighbors + REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: grid_max,grid_min + INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: nijk + INTEGER, DIMENSION(:,:,:,:), ALLOCATABLE :: ijk + +! *** Public subroutines *** + + PUBLIC :: build_fist_neighbor_lists,& + rebuild_fist_neighbor_lists + +! ***************************************************************************** + +CONTAINS + +! ***************************************************************************** + + SUBROUTINE build_fist_neighbor_lists(atomic_kind_set,particle_set,pnode,& + cell,r_cut,nonbonded,globenv) + +! Purpose: Build all the required neighbor lists for FIST. + +! History: - Creation (19.11.2002,MK) + +! *************************************************************************** + + TYPE(cell_type), POINTER :: cell + TYPE(global_environment_type), INTENT(IN) :: globenv + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: nonbonded + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(particle_node_type), DIMENSION(:), POINTER :: pnode + REAL(wp), DIMENSION(:,:), INTENT(IN) :: r_cut + +! *** Local parameters *** + + CHARACTER(LEN=*), PARAMETER :: routine_name = "build_fist_neighbor_lists" + CHARACTER(LEN=*), PARAMETER :: routine =& + "SUBROUTINE "//routine_name//" (MODULE "//module_name//")" + +! *** Local variables *** + + TYPE(atomic_kind_type), POINTER :: atomic_kind + INTEGER :: handle,nparticle + +! --------------------------------------------------------------------------- + + CALL write_checkpoint_information("entering "//routine_name,globenv) + + CALL timeset(routine_name,"I","",handle) + + group = globenv%group + ionode = globenv%ionode + mype = globenv%mepos + npe = globenv%num_pe + output_unit = globenv%scr + + IF ((ionode.AND.globenv%print%cell_parameters).AND.& + (globenv%print%level > LOW)) THEN + print_cell_parameters = .TRUE. + ELSE + print_cell_parameters = .FALSE. + END IF + + nparticle = SIZE(particle_set) + + ALLOCATE (atom_of_kind(nparticle),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"atom_of_kind",nparticle*int_size) + ALLOCATE (kind_of(nparticle),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"kind_of",nparticle*int_size) + + CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& + atom_of_kind=atom_of_kind,& + kind_of=kind_of,& + maxatom=maxatom) + + CALL get_cell(cell=cell,& + periodic=periodic,& + subcells=subcells,& + unit_of_length=unit_of_length,& + unit_of_length_name=unit_of_length_name) + +!MK *** pnode and neighbor list structure do not fit very well *** + + nnode = SIZE(pnode) + maxexcl = 0 + DO inode=1,nnode + maxexcl = MAX(maxexcl,pnode(inode)%nexcl) + END DO + +! *** Allocate work storage *** + + ALLOCATE (kind_a(maxatom),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"kind_a",maxatom*int_size) + + nkind = SIZE(atomic_kind_set) + + ALLOCATE (excl_a(nkind),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"excl_a",nkind*int_size) + + ALLOCATE (natom(nkind),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"natom",nkind*int_size) + + ALLOCATE (atoms(nkind),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"atoms",nkind*int_size) + + ALLOCATE (pbc_coord(nkind),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"pbc_coord",nkind*int_size) + +! *** Calculate PBC coordinates *** + + DO ikind=1,nkind + + atomic_kind => atomic_kind_set(ikind) + + NULLIFY (atoms(ikind)%list) + + CALL get_atomic_kind(atomic_kind=atomic_kind,& + natom=natom(ikind),& + atom_list=atoms(ikind)%list) + + ALLOCATE (excl_a(ikind)%list(maxexcl,natom(ikind)),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"excl_a(ikind)%list",& + maxexcl*natom(ikind)*wp_size) + excl_a(ikind)%list(:,:) = 0 + + ALLOCATE (pbc_coord(ikind)%r(3,natom(ikind)),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"pbc_coord(ikind)%r",& + 3*natom(ikind)*wp_size) + ALLOCATE (pbc_coord(ikind)%s(3,natom(ikind)),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"pbc_coord(ikind)%s",& + 3*natom(ikind)*wp_size) + + DO iatom=1,natom(ikind) + atom_a = atoms(ikind)%list(iatom) + ra_pbc(:) = pbc(particle_set(atom_a)%r(:),cell) + pbc_coord(ikind)%r(:,iatom) = ra_pbc(:) + pbc_coord(ikind)%s(:,iatom) = real_to_scaled(ra_pbc(:),cell) + END DO + + END DO + + DO inode=1,nnode + atom_a = pnode(inode)%p%iatom + ikind = kind_of(atom_a) + iatom = atom_of_kind(atom_a) + excl_node => pnode(inode)%ex + DO iexcl=1,pnode(inode)%nexcl + excl_a(ikind)%list(iexcl,iatom) = excl_node%p%iatom + excl_node => excl_node%next + END DO + END DO + +! *** Build the nonbonded neighbor lists *** + + CALL build_nonbonded(atomic_kind_set,particle_set,cell,r_cut,nonbonded) + + IF (ionode.AND.globenv%print%sab_orb_neighbor_lists) THEN + CALL write_neighbor_lists(nonbonded,"NONBONDED NEIGHBOR LISTS",& + particle_set,cell) + END IF + +! *** Release work storage *** + + DEALLOCATE (atom_of_kind,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"atom_of_kind") + + DEALLOCATE (kind_of,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"kind_of") + + DEALLOCATE (kind_a,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"kind_a") + + DEALLOCATE (natom,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"natom") + + DO ikind=1,nkind + NULLIFY (atoms(ikind)%list) + IF (ASSOCIATED(excl_a(ikind)%list)) THEN + DEALLOCATE (excl_a(ikind)%list,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"excl_a(ikind)%list") + END IF + IF (ASSOCIATED(pbc_coord(ikind)%r)) THEN + DEALLOCATE (pbc_coord(ikind)%r,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"pbc_coord(ikind)%r") + END IF + IF (ASSOCIATED(pbc_coord(ikind)%s)) THEN + DEALLOCATE (pbc_coord(ikind)%s,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"pbc_coord(ikind)%s") + END IF + END DO + + DEALLOCATE (atoms,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"atoms") + + DEALLOCATE (excl_a,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"excl_a") + + DEALLOCATE (pbc_coord,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"pbc_coord") + + CALL timestop(0.0_wp,handle) + + CALL write_checkpoint_information("leaving "//routine_name,globenv) + + END SUBROUTINE build_fist_neighbor_lists + +! ***************************************************************************** + + SUBROUTINE build_nonbonded(atomic_kind_set,particle_set,cell,r_cut,nonbonded) + +! Purpose: Build the nonbonded neighbor lists for FIST. + +! History: - Creation (19.11.2002,MK) + +! *************************************************************************** + + TYPE(cell_type), POINTER :: cell + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: nonbonded + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + REAL(wp), DIMENSION(:,:), INTENT(IN) :: r_cut + +! *** Local parameters *** + + CHARACTER(LEN=*), PARAMETER :: routine_name = "build_nonbonded" + CHARACTER(LEN=*), PARAMETER :: routine =& + "SUBROUTINE "//routine_name//" (MODULE "//module_name//")" + +! *** Local variables *** + + INTEGER :: ab,handle,natom_a + LOGICAL :: equal_kinds + +! --------------------------------------------------------------------------- + + IF (ASSOCIATED(nonbonded)) THEN + DO ab=1,SIZE(nonbonded) + CALL deallocate_neighbor_list_set(nonbonded(ab)%neighbor_list_set) + END DO + DEALLOCATE (nonbonded,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"nonbonded") + END IF + + ALLOCATE (nonbonded(nkind*(nkind + 1)/2),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"nonbonded",0) + + DO ab=1,SIZE(nonbonded) + NULLIFY (nonbonded(ab)%neighbor_list_set) + END DO + +! *** Print headline *** + + IF (print_cell_parameters) THEN + WRITE (UNIT=output_unit,FMT="(/,/,T2,A,/,/,T3,A,T29,A,T54,A)")& + "SUBCELL GRID FOR THE NONBONDED NEIGHBOR LISTS",& + "Atomic kind pair","Grid size",& + "Subcell size in "//unit_of_length_name + END IF + +! *** Loop over all atomic kind pairs *** + + DO ikind=1,nkind + DO jkind=ikind,nkind + + ab = ikind + jkind*(jkind - 1)/2 + + equal_kinds = (ikind == jkind) + +! *** Calculate the square of the maximum interaction distance *** + + rab_max = r_cut(ikind,jkind) + rab2_max = rab_max*rab_max + + r(:) = rab_max + sab_max(:) = real_to_scaled(r(:),cell) + + ncell(:) = (INT(sab_max(:)) + 1)*periodic(:) + ngrid(:) = MAX(1,NINT(0.5_wp*subcells/sab_max(:))) + +! *** Print subcell information for the current atomic kind pair *** + + IF (print_cell_parameters) THEN + WRITE (UNIT=output_unit,FMT="(T3,2I8,4X,3I5,6X,3F12.6)")& + ikind,jkind,ngrid,& + scaled_to_real(1.0_wp/REAL(ngrid(:),wp),cell)/unit_of_length + END IF + + CALL allocate_neighbor_list_set(neighbor_list_set=& + nonbonded(ab)%neighbor_list_set,& + r_max=rab_max) + neighbor_list_set => nonbonded(ab)%neighbor_list_set + + cell_a = (/0,0,0/) + +! *** Check, if we have to consider a subcell grid *** + + IF (SUM(ngrid) == 3) THEN + + DO iatom=1,natom(ikind) + atom_a = atoms(ikind)%list(iatom) + CALL add_neighbor_list(neighbor_list_set=neighbor_list_set,& + atom=atom_a,& + cell=cell_a,& + neighbor_list=kind_a(iatom)%neighbor_list) + END DO + + DO jatom=1,natom(jkind) + + atom_b = atoms(jkind)%list(jatom) + sb_pbc(:) = pbc_coord(jkind)%s(:,jatom) + + DO icell=-ncell(1),ncell(1) + cell_b(1) = icell + DO jcell=-ncell(2),ncell(2) + cell_b(2) = jcell + DO kcell=-ncell(3),ncell(3) + cell_b(3) = kcell + + sb(:) = sb_pbc(:) + REAL(cell_b(:),wp) + rb(:) = scaled_to_real(sb(:),cell) + + IF (equal_kinds) THEN + natom_a = jatom + ELSE + natom_a = natom(ikind) + END IF + + DO iatom=1,natom_a + rab(:) = rb(:) - pbc_coord(ikind)%r(:,iatom) + rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) + IF (rab2 < rab2_max) THEN + IF (rab2 > 1.0E-6_wp) THEN + CALL add_neighbor_node(& + neighbor_list=kind_a(iatom)%neighbor_list,& + neighbor=atom_b,& + cell=cell_b,& + r=rab(:),& + exclusion_list=excl_a(ikind)%list(:,iatom)) + END IF + END IF + END DO + + END DO + END DO + END DO + + END DO + + ELSE + + ALLOCATE (grid_min(3,ngrid(1),ngrid(2),ngrid(3)),STAT=istat) + IF (istat /= 0) THEN + CALL stop_memory(routine,"grid_min",3*PRODUCT(ngrid)*wp_size) + END IF + + ALLOCATE (grid_max(3,ngrid(1),ngrid(2),ngrid(3)),STAT=istat) + IF (istat /= 0) THEN + CALL stop_memory(routine,"grid_max",3*PRODUCT(ngrid)*wp_size) + END IF + + ALLOCATE (nijk(ngrid(1),ngrid(2),ngrid(3)),STAT=istat) + IF (istat /= 0) THEN + CALL stop_memory(routine,"nijk",PRODUCT(ngrid)*int_size) + END IF + nijk(:,:,:) = 0 + + ALLOCATE (ijk(natom(ikind),ngrid(1),ngrid(2),ngrid(3)),STAT=istat) + IF (istat /= 0) THEN + CALL stop_memory(routine,"ijk",& + natom(ikind)*PRODUCT(ngrid)*int_size) + END IF + + DO igrid=1,ngrid(1) + a_min = REAL(igrid-1,wp)/REAL(ngrid(1),wp) - 0.5_wp + a_max = REAL(igrid,wp)/REAL(ngrid(1),wp) - 0.5_wp + DO jgrid=1,ngrid(2) + b_min = REAL(jgrid-1,wp)/REAL(ngrid(2),wp) - 0.5_wp + b_max = REAL(jgrid,wp)/REAL(ngrid(2),wp) - 0.5_wp + DO kgrid=1,ngrid(3) + c_min = REAL(kgrid-1,wp)/REAL(ngrid(3),wp) - 0.5_wp + c_max = REAL(kgrid,wp)/REAL(ngrid(3),wp) - 0.5_wp + grid_min(:,igrid,jgrid,kgrid) = (/a_min,b_min,c_min/) + grid_max(:,igrid,jgrid,kgrid) = (/a_max,b_max,c_max/) + END DO + END DO + END DO + + DO iatom=1,natom(ikind) + atom_a = atoms(ikind)%list(iatom) + sa_pbc(:) = pbc_coord(ikind)%s(:,iatom) + igrid = MAX(1,CEILING((sa_pbc(1) + 0.5_wp)*ngrid(1))) + jgrid = MAX(1,CEILING((sa_pbc(2) + 0.5_wp)*ngrid(2))) + kgrid = MAX(1,CEILING((sa_pbc(3) + 0.5_wp)*ngrid(3))) + nijk(igrid,jgrid,kgrid) = nijk(igrid,jgrid,kgrid) + 1 + ijk(nijk(igrid,jgrid,kgrid),igrid,jgrid,kgrid) = iatom + CALL add_neighbor_list(neighbor_list_set=neighbor_list_set,& + atom=atom_a,& + cell=cell_a,& + neighbor_list=kind_a(iatom)%neighbor_list) + END DO + + DO jatom=1,natom(jkind) + + atom_b = atoms(jkind)%list(jatom) + sb_pbc(:) = pbc_coord(jkind)%s(:,jatom) + + DO icell=-ncell(1),ncell(1) + cell_b(1) = icell + DO jcell=-ncell(2),ncell(2) + cell_b(2) = jcell + DO kcell=-ncell(3),ncell(3) + cell_b(3) = kcell + + sb(:) = sb_pbc(:) + REAL(cell_b(:),wp) + rb(:) = scaled_to_real(sb(:),cell) + sb_min(:) = sb(:) - sab_max(:) + sb_max(:) = sb(:) + sab_max(:) + + IF (sb_max(1) < grid_min(1,1,1,1)) CYCLE + IF (sb_max(2) < grid_min(2,1,1,1)) CYCLE + IF (sb_max(3) < grid_min(3,1,1,1)) CYCLE + + igrid = ngrid(1) + jgrid = ngrid(2) + kgrid = ngrid(3) + + IF (sb_min(1) >= grid_max(1,igrid,jgrid,kgrid)) CYCLE + IF (sb_min(2) >= grid_max(2,igrid,jgrid,kgrid)) CYCLE + IF (sb_min(3) >= grid_max(3,igrid,jgrid,kgrid)) CYCLE + + DO igrid=1,ngrid(1) + DO jgrid=1,ngrid(2) + DO kgrid=1,ngrid(3) + + IF (nijk(igrid,jgrid,kgrid) == 0) CYCLE + + IF (sb_max(1) < grid_min(1,igrid,jgrid,kgrid)) CYCLE + IF (sb_max(2) < grid_min(2,igrid,jgrid,kgrid)) CYCLE + IF (sb_max(3) < grid_min(3,igrid,jgrid,kgrid)) CYCLE + + IF (sb_min(1) >= grid_max(1,igrid,jgrid,kgrid)) CYCLE + IF (sb_min(2) >= grid_max(2,igrid,jgrid,kgrid)) CYCLE + IF (sb_min(3) >= grid_max(3,igrid,jgrid,kgrid)) CYCLE + + DO iijk=1,nijk(igrid,jgrid,kgrid) + iatom = ijk(iijk,igrid,jgrid,kgrid) + IF (equal_kinds) THEN + IF (jatom < iatom) CYCLE + END IF + rab(:) = rb(:) - pbc_coord(ikind)%r(:,iatom) + rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) + IF (rab2 < rab2_max) THEN + IF (rab2 > 1.0E-6_wp) THEN + CALL add_neighbor_node(& + neighbor_list=kind_a(iatom)%neighbor_list,& + neighbor=atom_b,& + cell=cell_b,& + r=rab(:),& + exclusion_list=excl_a(ikind)%list(:,iatom)) + END IF + END IF + END DO + + END DO + END DO + END DO + + END DO + END DO + END DO + + END DO + + DEALLOCATE (grid_min,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"grid_min") + + DEALLOCATE (grid_max,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"grid_max") + + DEALLOCATE (nijk,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"nijk") + + DEALLOCATE (ijk,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"ijk") + + END IF + + END DO + END DO + + END SUBROUTINE build_nonbonded + +! ***************************************************************************** + + SUBROUTINE rebuild_fist_neighbor_lists(atomic_kind_set,particle_set,pnode,& + cell,r_cut,nonbonded,globenv) + +! Purpose: Rebuild all the required neighbor lists for FIST. + +! History: - Creation (19.11.2002,MK) + +! *************************************************************************** + + TYPE(cell_type), POINTER :: cell + TYPE(global_environment_type), INTENT(IN) :: globenv + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: nonbonded + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + TYPE(particle_node_type), DIMENSION(:), POINTER :: pnode + REAL(wp), DIMENSION(:,:), INTENT(IN) :: r_cut + +! *** Local parameters *** + + CHARACTER(LEN=*), PARAMETER :: routine_name = "rebuild_fist_neighbor_lists" + CHARACTER(LEN=*), PARAMETER :: routine =& + "SUBROUTINE "//routine_name//" (MODULE "//module_name//")" + +! *** Local variables *** + + TYPE(atomic_kind_type), POINTER :: atomic_kind + + INTEGER :: handle,nparticle + +! --------------------------------------------------------------------------- + + CALL write_checkpoint_information("entering "//routine_name,globenv) + + CALL timeset(routine_name,"I","",handle) + + group = globenv%group + ionode = globenv%ionode + mype = globenv%mepos + npe = globenv%num_pe + output_unit = globenv%scr + + nparticle = SIZE(particle_set) + + ALLOCATE (atom_of_kind(nparticle),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"atom_of_kind",nparticle*int_size) + ALLOCATE (kind_of(nparticle),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"kind_of",nparticle*int_size) + + CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set,& + atom_of_kind=atom_of_kind,& + kind_of=kind_of,& + maxatom=maxatom) + + CALL get_cell(cell=cell,& + periodic=periodic,& + subcells=subcells,& + unit_of_length=unit_of_length,& + unit_of_length_name=unit_of_length_name) + +!MK *** pnode and neighbor list structure do not fit very well *** + + nnode = SIZE(pnode) + maxexcl = 0 + DO inode=1,nnode + maxexcl = MAX(maxexcl,pnode(inode)%nexcl) + END DO + +! *** Allocate work storage *** + + ALLOCATE (kind_a(maxatom),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"kind_a",maxatom*int_size) + + nkind = SIZE(atomic_kind_set) + + ALLOCATE (excl_a(nkind),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"excl_a",nkind*int_size) + + ALLOCATE (natom(nkind),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"natom",nkind*int_size) + + ALLOCATE (atoms(nkind),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"atoms",nkind*int_size) + + ALLOCATE (pbc_coord(nkind),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"pbc_coord",nkind*int_size) + +! *** Calculate PBC coordinates *** + + DO ikind=1,nkind + + atomic_kind => atomic_kind_set(ikind) + + NULLIFY (atoms(ikind)%list) + + CALL get_atomic_kind(atomic_kind=atomic_kind,& + natom=natom(ikind),& + atom_list=atoms(ikind)%list) + + ALLOCATE (excl_a(ikind)%list(maxexcl,natom(ikind)),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"excl_a(ikind)%list",& + maxexcl*natom(ikind)*wp_size) + excl_a(ikind)%list(:,:) = 0 + + ALLOCATE (pbc_coord(ikind)%r(3,natom(ikind)),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"pbc_coord(ikind)%r",& + 3*natom(ikind)*wp_size) + ALLOCATE (pbc_coord(ikind)%s(3,natom(ikind)),STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"pbc_coord(ikind)%s",& + 3*natom(ikind)*wp_size) + + DO iatom=1,natom(ikind) + atom_a = atoms(ikind)%list(iatom) + ra_pbc(:) = pbc(particle_set(atom_a)%r(:),cell) + pbc_coord(ikind)%r(:,iatom) = ra_pbc(:) + pbc_coord(ikind)%s(:,iatom) = real_to_scaled(ra_pbc(:),cell) + END DO + + END DO + + DO inode=1,nnode + atom_a = pnode(inode)%p%iatom + ikind = kind_of(atom_a) + iatom = atom_of_kind(atom_a) + excl_node => pnode(inode)%ex + DO iexcl=1,pnode(inode)%nexcl + excl_a(ikind)%list(iexcl,iatom) = excl_node%p%iatom + excl_node => excl_node%next + END DO + END DO + +! *** Rebuild the nonbonded neighbor lists *** + + CALL rebuild_nonbonded(atomic_kind_set,particle_set,cell,r_cut,nonbonded) + + IF (ionode.AND.globenv%print%sab_orb_neighbor_lists) THEN + CALL write_neighbor_lists(nonbonded,"NONBONDED NEIGHBOR LISTS",& + particle_set,cell) + END IF + +! *** Release work storage *** + + DEALLOCATE (atom_of_kind,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"atom_of_kind") + + DEALLOCATE (kind_of,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"kind_of") + + DEALLOCATE (kind_a,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"kind_a") + + DEALLOCATE (natom,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"natom") + + DO ikind=1,nkind + NULLIFY (atoms(ikind)%list) + IF (ASSOCIATED(excl_a(ikind)%list)) THEN + DEALLOCATE (excl_a(ikind)%list,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"excl_a(ikind)%list") + END IF + IF (ASSOCIATED(pbc_coord(ikind)%r)) THEN + DEALLOCATE (pbc_coord(ikind)%r,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"pbc_coord(ikind)%r") + END IF + IF (ASSOCIATED(pbc_coord(ikind)%s)) THEN + DEALLOCATE (pbc_coord(ikind)%s,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"pbc_coord(ikind)%s") + END IF + END DO + + DEALLOCATE (atoms,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"atoms") + + DEALLOCATE (excl_a,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"excl_a") + + DEALLOCATE (pbc_coord,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"pbc_coord") + + CALL timestop(0.0_wp,handle) + + CALL write_checkpoint_information("leaving "//routine_name,globenv) + + END SUBROUTINE rebuild_fist_neighbor_lists + +! ***************************************************************************** + + SUBROUTINE rebuild_nonbonded(atomic_kind_set,particle_set,cell,r_cut,& + nonbonded) + +! Purpose: Rebuild the nonbonded neighbor lists for FIST. + +! History: - Creation (19.11.2002,MK) + +! *************************************************************************** + + TYPE(cell_type), POINTER :: cell + TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set + TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: nonbonded + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + REAL(wp), DIMENSION(:,:), INTENT(IN) :: r_cut + +! *** Local parameters *** + + CHARACTER(LEN=*), PARAMETER :: routine_name = "rebuild_nonbonded" + CHARACTER(LEN=*), PARAMETER :: routine =& + "SUBROUTINE "//routine_name//" (MODULE "//module_name//")" + +! *** Local variables *** + + INTEGER :: ab,handle,natom_a + LOGICAL :: equal_kinds + +! --------------------------------------------------------------------------- + +! *** Loop over all atomic kind pairs *** + + DO ikind=1,nkind + DO jkind=ikind,nkind + + ab = ikind + jkind*(jkind - 1)/2 + + neighbor_list_set => nonbonded(ab)%neighbor_list_set + + IF (.NOT.ASSOCIATED(neighbor_list_set)) CYCLE + + equal_kinds = (ikind == jkind) + + rab_max = r_cut(ikind,jkind) + rab2_max = rab_max*rab_max + + r(:) = rab_max + sab_max(:) = real_to_scaled(r(:),cell) + + ncell(:) = (INT(sab_max(:)) + 1)*periodic(:) + ngrid(:) = MAX(1,NINT(0.5_wp*subcells/sab_max(:))) + + cell_a = (/0,0,0/) + +! *** Check, if we have to consider a subcell grid *** + + IF (SUM(ngrid) == 3) THEN + + neighbor_list => first_list(neighbor_list_set) + + DO WHILE (ASSOCIATED(neighbor_list)) + CALL get_neighbor_list(neighbor_list=neighbor_list,& + atom=atom_a) + CALL init_neighbor_list(neighbor_list) + iatom = atom_of_kind(atom_a) + kind_a(iatom)%neighbor_list => neighbor_list + neighbor_list => next(neighbor_list) + END DO + + DO jatom=1,natom(jkind) + + atom_b = atoms(jkind)%list(jatom) + sb_pbc(:) = pbc_coord(jkind)%s(:,jatom) + + DO icell=-ncell(1),ncell(1) + cell_b(1) = icell + DO jcell=-ncell(2),ncell(2) + cell_b(2) = jcell + DO kcell=-ncell(3),ncell(3) + cell_b(3) = kcell + + sb(:) = sb_pbc(:) + REAL(cell_b(:),wp) + rb(:) = scaled_to_real(sb(:),cell) + + IF (equal_kinds) THEN + natom_a = jatom + ELSE + natom_a = natom(ikind) + END IF + + DO iatom=1,natom_a + rab(:) = rb(:) - pbc_coord(ikind)%r(:,iatom) + rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) + IF (rab2 < rab2_max) THEN + IF (rab2 > 1.0E-6_wp) THEN + CALL add_neighbor_node(& + neighbor_list=kind_a(iatom)%neighbor_list,& + neighbor=atom_b,& + cell=cell_b,& + r=rab(:),& + exclusion_list=excl_a(ikind)%list(:,iatom)) + END IF + END IF + END DO + + END DO + END DO + END DO + + END DO + + ELSE + + ALLOCATE (grid_min(3,ngrid(1),ngrid(2),ngrid(3)),STAT=istat) + IF (istat /= 0) THEN + CALL stop_memory(routine,"grid_min",3*PRODUCT(ngrid)*wp_size) + END IF + + ALLOCATE (grid_max(3,ngrid(1),ngrid(2),ngrid(3)),STAT=istat) + IF (istat /= 0) THEN + CALL stop_memory(routine,"grid_max",3*PRODUCT(ngrid)*wp_size) + END IF + + ALLOCATE (nijk(ngrid(1),ngrid(2),ngrid(3)),STAT=istat) + IF (istat /= 0) THEN + CALL stop_memory(routine,"nijk",PRODUCT(ngrid)*int_size) + END IF + nijk(:,:,:) = 0 + + ALLOCATE (ijk(natom(ikind),ngrid(1),ngrid(2),ngrid(3)),STAT=istat) + IF (istat /= 0) THEN + CALL stop_memory(routine,"ijk",& + natom(ikind)*PRODUCT(ngrid)*int_size) + END IF + + DO igrid=1,ngrid(1) + a_min = REAL(igrid-1,wp)/REAL(ngrid(1),wp) - 0.5_wp + a_max = REAL(igrid,wp)/REAL(ngrid(1),wp) - 0.5_wp + DO jgrid=1,ngrid(2) + b_min = REAL(jgrid-1,wp)/REAL(ngrid(2),wp) - 0.5_wp + b_max = REAL(jgrid,wp)/REAL(ngrid(2),wp) - 0.5_wp + DO kgrid=1,ngrid(3) + c_min = REAL(kgrid-1,wp)/REAL(ngrid(3),wp) - 0.5_wp + c_max = REAL(kgrid,wp)/REAL(ngrid(3),wp) - 0.5_wp + grid_min(:,igrid,jgrid,kgrid) = (/a_min,b_min,c_min/) + grid_max(:,igrid,jgrid,kgrid) = (/a_max,b_max,c_max/) + END DO + END DO + END DO + + neighbor_list => first_list(neighbor_list_set) + + DO WHILE (ASSOCIATED(neighbor_list)) + CALL get_neighbor_list(neighbor_list=neighbor_list,& + atom=atom_a) + CALL init_neighbor_list(neighbor_list) + iatom = atom_of_kind(atom_a) + sa_pbc(:) = pbc_coord(ikind)%s(:,iatom) + igrid = MAX(1,CEILING((sa_pbc(1) + 0.5_wp)*ngrid(1))) + jgrid = MAX(1,CEILING((sa_pbc(2) + 0.5_wp)*ngrid(2))) + kgrid = MAX(1,CEILING((sa_pbc(3) + 0.5_wp)*ngrid(3))) + nijk(igrid,jgrid,kgrid) = nijk(igrid,jgrid,kgrid) + 1 + ijk(nijk(igrid,jgrid,kgrid),igrid,jgrid,kgrid) = iatom + kind_a(iatom)%neighbor_list => neighbor_list + neighbor_list => next(neighbor_list) + END DO + + DO jatom=1,natom(jkind) + + atom_b = atoms(jkind)%list(jatom) + sb_pbc(:) = pbc_coord(jkind)%s(:,jatom) + + DO icell=-ncell(1),ncell(1) + cell_b(1) = icell + DO jcell=-ncell(2),ncell(2) + cell_b(2) = jcell + DO kcell=-ncell(3),ncell(3) + cell_b(3) = kcell + + sb(:) = sb_pbc(:) + REAL(cell_b(:),wp) + rb(:) = scaled_to_real(sb(:),cell) + sb_min(:) = sb(:) - sab_max(:) + sb_max(:) = sb(:) + sab_max(:) + + IF (sb_max(1) < grid_min(1,1,1,1)) CYCLE + IF (sb_max(2) < grid_min(2,1,1,1)) CYCLE + IF (sb_max(3) < grid_min(3,1,1,1)) CYCLE + + igrid = ngrid(1) + jgrid = ngrid(2) + kgrid = ngrid(3) + + IF (sb_min(1) >= grid_max(1,igrid,jgrid,kgrid)) CYCLE + IF (sb_min(2) >= grid_max(2,igrid,jgrid,kgrid)) CYCLE + IF (sb_min(3) >= grid_max(3,igrid,jgrid,kgrid)) CYCLE + + DO igrid=1,ngrid(1) + DO jgrid=1,ngrid(2) + DO kgrid=1,ngrid(3) + + IF (nijk(igrid,jgrid,kgrid) == 0) CYCLE + + IF (sb_max(1) < grid_min(1,igrid,jgrid,kgrid)) CYCLE + IF (sb_max(2) < grid_min(2,igrid,jgrid,kgrid)) CYCLE + IF (sb_max(3) < grid_min(3,igrid,jgrid,kgrid)) CYCLE + + IF (sb_min(1) >= grid_max(1,igrid,jgrid,kgrid)) CYCLE + IF (sb_min(2) >= grid_max(2,igrid,jgrid,kgrid)) CYCLE + IF (sb_min(3) >= grid_max(3,igrid,jgrid,kgrid)) CYCLE + + DO iijk=1,nijk(igrid,jgrid,kgrid) + iatom = ijk(iijk,igrid,jgrid,kgrid) + IF (equal_kinds) THEN + IF (jatom < iatom) CYCLE + END IF + rab(:) = rb(:) - pbc_coord(ikind)%r(:,iatom) + rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3) + IF (rab2 < rab2_max) THEN + IF (rab2 > 1.0E-6_wp) THEN + CALL add_neighbor_node(& + neighbor_list=kind_a(iatom)%neighbor_list,& + neighbor=atom_b,& + cell=cell_b,& + r=rab(:),& + exclusion_list=excl_a(ikind)%list(:,iatom)) + END IF + END IF + END DO + + END DO + END DO + END DO + + END DO + END DO + END DO + + END DO + +! *** Release work storage *** + + DEALLOCATE (grid_min,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"grid_min") + + DEALLOCATE (grid_max,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"grid_max") + + DEALLOCATE (nijk,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"nijk") + + DEALLOCATE (ijk,STAT=istat) + IF (istat /= 0) CALL stop_memory(routine,"ijk") + + END IF + +!MK CALL clean_neighbor_list_set(neighbor_list_set) + + END DO + END DO + + END SUBROUTINE rebuild_nonbonded + +! ***************************************************************************** + + SUBROUTINE write_neighbor_lists(neighbor_lists,name,particle_set,cell) + +! Purpose: Write a set of neighbor lists to the output unit. + +! History: - Creation (04.03.2002,MK) + +! *************************************************************************** + + TYPE(cell_type), POINTER :: cell + TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: neighbor_lists + TYPE(particle_type), DIMENSION(:), POINTER :: particle_set + CHARACTER(LEN=*), INTENT(IN) :: name + +! *** Local variables *** + + INTEGER :: ab,atom_a,atom_b,i,nneighbor + + REAL(wp), DIMENSION(3) :: ra,rab,rb + INTEGER, DIMENSION(3) :: cell_a,cell_b + +! --------------------------------------------------------------------------- + +! *** Write headline *** + + WRITE (UNIT=output_unit,FMT="(/,/,T2,A,/,/,T3,A,7X,A,2(11X,A),10X,A)")& + TRIM(name)//" IN "//TRIM(unit_of_length_name),& + "Atom Neighbor Cell(i,j,k)","X","Y","Z","Distance" + + DO ab=1,SIZE(neighbor_lists) + + neighbor_list_set => neighbor_lists(ab)%neighbor_list_set + + IF (.NOT.ASSOCIATED(neighbor_list_set)) CYCLE + +! *** Loop over all atoms and their corresponding neighbor lists *** + + neighbor_list => first_list(neighbor_list_set) + + DO WHILE (ASSOCIATED(neighbor_list)) + + CALL get_neighbor_list(neighbor_list=neighbor_list,& + atom=atom_a,& + cell=cell_a,& + nnode=nneighbor) + + ra(:) = pbc(particle_set(atom_a)%r,cell,cell_a) + + WRITE (UNIT=output_unit,FMT="(/,T2,I5,3X,I6,2X,3I4,3F12.6)")& + atom_a,nneighbor,cell_a(:),ra(:)/unit_of_length + +! *** Direct the work pointer to the *** +! *** start point of the current list *** + + neighbor_node => first_node(neighbor_list) + +! *** Traverse the neighbor list of the current *** +! *** atom and print the stored information *** + + DO WHILE (ASSOCIATED(neighbor_node)) + + CALL get_neighbor_node(neighbor_node=neighbor_node,& + neighbor=atom_b,& + cell=cell_b,& + r=rab) + + rb(:) = ra(:) + rab(:) + + WRITE (UNIT=output_unit,FMT="(T10,I6,2X,3I4,3F12.6,2X,F12.6)")& + atom_b,cell_b,rb(:)/unit_of_length,& + SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))/unit_of_length + + neighbor_node => next(neighbor_node) + + END DO + + neighbor_list => next(neighbor_list) + + END DO + + END DO + + END SUBROUTINE write_neighbor_lists + +! ***************************************************************************** + +END MODULE fist_neighbor_lists diff --git a/src/fist_nonbond_force.F b/src/fist_nonbond_force.F index c8885be901..067bffae5a 100644 --- a/src/fist_nonbond_force.F +++ b/src/fist_nonbond_force.F @@ -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 diff --git a/src/linklist_control.F b/src/linklist_control.F index acddf0620d..7888eacf55 100644 --- a/src/linklist_control.F +++ b/src/linklist_control.F @@ -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 diff --git a/src/linklist_types.F b/src/linklist_types.F index fb273080ed..fcff0ac5c2 100644 --- a/src/linklist_types.F +++ b/src/linklist_types.F @@ -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 diff --git a/src/linklist_verlet_list.F b/src/linklist_verlet_list.F index d3495a04ce..c5f37ea9b3 100644 --- a/src/linklist_verlet_list.F +++ b/src/linklist_verlet_list.F @@ -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 diff --git a/src/md.F b/src/md.F index 1d4ee01435..985f458efa 100644 --- a/src/md.F +++ b/src/md.F @@ -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 diff --git a/src/qs_neighbor_list_types.F b/src/qs_neighbor_list_types.F index a350cd8cb4..d05234eb99 100644 --- a/src/qs_neighbor_list_types.F +++ b/src/qs_neighbor_list_types.F @@ -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)