diff --git a/AUTHORS b/AUTHORS index 59bd0d1..4e9ad04 100644 --- a/AUTHORS +++ b/AUTHORS @@ -9,3 +9,8 @@ dynamics and docking. Journal of Computational Chemistry, 24(7):819-825, 2003. Home-Page: http://m3g.iqm.unicamp.br/packmol + +Contributor: +- Michele Bonus, Heinrich Heine University Düsseldorf + Email: michele.bonus@hhu.de + GitHub: https://github.com/MicheleBonus/ diff --git a/CHANGELOG.md b/CHANGELOG.md index efbc3bd..48ca224 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,8 @@ packmol changelog Version 21.2.2-DEV -------------- - ![INFO][badge-info] Add changelog verification CI run. +- ![INFO][badge-info] Update README build/profile instructions and document recent cell-list and pair-kernel performance work. +- ![INFO][badge-info] Add contributor credentials for Michele Bonus (Heinrich Heine University Düsseldorf). Version 21.2.1 -------------- diff --git a/CMakeLists.txt b/CMakeLists.txt index 860b74b..149dc3b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,15 +1,48 @@ cmake_minimum_required(VERSION 3.10) project(packmol Fortran) +set(PACKMOL_PROFILE "baseline" CACHE STRING "Build profile: baseline, perf-native, debug, sanitize") +set_property(CACHE PACKMOL_PROFILE PROPERTY STRINGS baseline perf-native debug sanitize) +option(PACKMOL_UNSAFE_MATH "Enable unsafe/fast math optimizations" OFF) if(NOT CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE Release) -endif(NOT CMAKE_BUILD_TYPE) +endif() if(CMAKE_Fortran_COMPILER_ID MATCHES GNU) - add_compile_options( - -Wall "$<$:-Werror>" - ) + add_compile_options(-Wall) +endif() + +if(PACKMOL_PROFILE STREQUAL "baseline") + add_compile_options(-O2) +elseif(PACKMOL_PROFILE STREQUAL "perf-native") + add_compile_options(-O3) + if(CMAKE_Fortran_COMPILER_ID MATCHES GNU) + add_compile_options(-march=native -mtune=native -funroll-loops) + elseif(CMAKE_Fortran_COMPILER_ID MATCHES IntelLLVM OR CMAKE_Fortran_COMPILER_ID MATCHES Intel) + add_compile_options(-xHost -ipo) + endif() +elseif(PACKMOL_PROFILE STREQUAL "debug") + add_compile_options(-g) + if(CMAKE_Fortran_COMPILER_ID MATCHES GNU) + add_compile_options(-fcheck=bounds -fbacktrace -ffpe-trap=zero,overflow,underflow) + endif() +elseif(PACKMOL_PROFILE STREQUAL "sanitize") + add_compile_options(-O1 -g) + if(CMAKE_Fortran_COMPILER_ID MATCHES GNU) + add_compile_options(-fno-omit-frame-pointer -fsanitize=address,undefined -fcheck=all) + add_link_options(-fsanitize=address,undefined) + endif() +else() + message(FATAL_ERROR "Unknown PACKMOL_PROFILE='${PACKMOL_PROFILE}'. Use one of: baseline, perf-native, debug, sanitize") +endif() + +if(PACKMOL_UNSAFE_MATH) + if(CMAKE_Fortran_COMPILER_ID MATCHES GNU) + add_compile_options(-ffast-math) + elseif(CMAKE_Fortran_COMPILER_ID MATCHES IntelLLVM OR CMAKE_Fortran_COMPILER_ID MATCHES Intel) + add_compile_options(-fp-model fast=2) + endif() endif() # Build the executable diff --git a/Makefile b/Makefile index 11a0535..b0db730 100644 --- a/Makefile +++ b/Makefile @@ -18,7 +18,26 @@ FORTRAN=/usr/bin/gfortran # # Change the flags of the compilation if you want: # -FLAGS= -O3 -march=native -funroll-loops +PROFILE ?= baseline +UNSAFE_MATH ?= 0 + +ifneq (,$(filter perf perf-native,$(MAKECMDGOALS))) +PROFILE := perf-native +endif +ifneq (,$(filter devel,$(MAKECMDGOALS))) +PROFILE := devel +endif +ifneq (,$(filter sanitize,$(MAKECMDGOALS))) +PROFILE := sanitize +endif +ifneq (,$(filter static,$(MAKECMDGOALS))) +PROFILE := static +endif +ifneq (,$(filter baseline,$(MAKECMDGOALS))) +PROFILE := baseline +endif + +FLAGS= SRCDIR= src MAINDIR= app ################################################################### @@ -30,20 +49,54 @@ MAINDIR= app # Flags for compiling development version # GENCANFLAGS := $(FLAGS) -# Flags for the routines that signal with --fast-math -IEEE_SIGNAL_FLAGS := $(FLAGS) -ifeq ($(MAKECMDGOALS),devel) +FC_VERSION := $(shell $(FORTRAN) --version 2>/dev/null | head -n 1) + +ifeq ($(PROFILE),baseline) +FLAGS = -O2 +GENCANFLAGS = $(FLAGS) +endif + +ifeq ($(PROFILE),perf-native) +FLAGS = -O3 +GENCANFLAGS = $(FLAGS) +ifneq (,$(findstring GNU Fortran,$(FC_VERSION))) +FLAGS += -march=native -mtune=native -funroll-loops +GENCANFLAGS += -march=native -mtune=native -funroll-loops +else ifneq (,$(findstring ifx,$(FC_VERSION))) +FLAGS += -xHost -ipo +GENCANFLAGS += -xHost -ipo +else ifneq (,$(findstring ifort,$(FC_VERSION))) +FLAGS += -xHost -ipo +GENCANFLAGS += -xHost -ipo +endif +endif + +ifeq ($(PROFILE),devel) FLAGS = -Wall -fcheck=bounds -g -fbacktrace -ffpe-trap=zero,overflow,underflow -GENCANFLAGS = -fcheck=bounds -g -fbacktrace -ffpe-trap=zero,overflow,underflow +GENCANFLAGS = -fcheck=bounds -g -fbacktrace -ffpe-trap=zero,overflow,underflow endif -ifeq ($(MAKECMDGOALS),perf) -FLAGS = -g -pg -GENCANFLAGS = -g -pg + +ifeq ($(PROFILE),sanitize) +FLAGS = -O1 -g -fno-omit-frame-pointer +GENCANFLAGS = $(FLAGS) +ifneq (,$(findstring GNU Fortran,$(FC_VERSION))) +FLAGS += -fsanitize=address,undefined -fcheck=all +GENCANFLAGS += -fsanitize=address,undefined -fcheck=all endif -ifeq ($(MAKECMDGOALS),static) -FLAGS = -O3 --fast-math -static -GENCANFLAGS = -O3 --fast-math -static endif + +ifeq ($(PROFILE),static) +FLAGS = -O2 -static +GENCANFLAGS = $(FLAGS) +endif + +ifeq ($(UNSAFE_MATH),1) +FLAGS += -ffast-math +GENCANFLAGS += -ffast-math +endif +IEEE_SIGNAL_FLAGS = $(FLAGS) + + # # Files required # @@ -103,17 +156,22 @@ all : $(oall) # # Compiling with flags for development # -static : devel -perf : devel -devel : $(oall) - @echo " ------------------------------------------------------ " - @echo " Compiling packmol with $(FORTRAN) " - @echo " Flags: $(FLAGS)" - @echo " ------------------------------------------------------ " - @$(FORTRAN) -o packmol $(oall) $(FLAGS) - @echo " ------------------------------------------------------ " - @echo " Packmol succesfully built. " - @echo " ------------------------------------------------------ " +baseline: PROFILE=baseline +baseline: all + +perf-native: PROFILE=perf-native +perf-native: all + +devel: PROFILE=devel +devel: all + +sanitize: PROFILE=sanitize +sanitize: all + +static: PROFILE=static +static: all + +perf : perf-native # # Modules # diff --git a/README.md b/README.md index 79ef999..648ffcb 100644 --- a/README.md +++ b/README.md @@ -83,6 +83,101 @@ uvx packmol < inp.pack `FPM_FC=compiler`, for example for `ifort`, use in bash, `export FPM_FC=ifort`. + +### Build profiles (Makefile and CMake) + +Packmol now supports explicit build profiles so you can choose between conservative optimized builds, aggressive benchmarking builds, and debugging/sanitized builds. + +#### Baseline release profile (default) + +This profile prioritizes stability and portability with conservative optimization. + +- **Makefile**: `make` (or `make baseline`) +- **CMake**: `cmake -S . -B build -DPACKMOL_PROFILE=baseline` + +Then build with: + +```bash +cmake --build build --config Release +``` + +#### Aggressive benchmark profile (`perf-native`) + +This profile is opt-in and applies compiler-specific native tuning (for example, `-march=native` for GNU or `-xHost` for Intel compilers). Use this for local benchmarking only. + +- **Makefile**: `make perf-native` +- **CMake**: `cmake -S . -B build -DPACKMOL_PROFILE=perf-native` + +Then build with: + +```bash +cmake --build build --config Release +``` + +#### Debug / sanitized profiles + +Use these profiles for troubleshooting and memory/UB diagnostics. + +- **Makefile debug checks**: `make devel` +- **Makefile sanitizers (GNU)**: `make sanitize` +- **CMake debug checks**: `cmake -S . -B build -DPACKMOL_PROFILE=debug` +- **CMake sanitizers (GNU)**: `cmake -S . -B build -DPACKMOL_PROFILE=sanitize` + +Then build with: + +```bash +cmake --build build --config Debug +``` + +### Performance-oriented kernel improvements + +Recent performance work in the pairwise cell-list kernels includes: + +- Prebuilt neighbor-cell offsets used in `computef` / `computeg` traversal. +- Reduced branch pressure in pair kernels by separating hot paths. +- Better data locality in hot loops through structure-of-arrays style access. +- Conservative per-cell broad-phase reject checks before expensive pair evaluation. +- Consistent `x*x` math micro-optimizations replacing hot-loop `**2` exponentiation. +- Build-profile driven tuning with numerics validation support. + +These updates are designed to speed up large systems while preserving numerical behavior. + +## Contributors + +- Michele Bonus, Heinrich Heine University Düsseldorf + Email: `michele.bonus@hhu.de` + GitHub: https://github.com/MicheleBonus/ + +#### Unsafe math is explicit opt-in + +Unsafe/fast-math is **not** enabled globally. + +- **Makefile**: add `UNSAFE_MATH=1`, e.g. `make perf-native UNSAFE_MATH=1` +- **CMake**: set `-DPACKMOL_UNSAFE_MATH=ON` + +### Reproducibility and numerics check across profiles + +To compare objective and gradient behavior across profiles using existing inputs under `testing/`, run: + +```bash +./testing/check_numerics_profiles.sh +``` + +What it does: + +1. Builds and runs `baseline` and `perf-native` binaries. +2. Uses `testing/input_files/benzene2.inp` with `chkgrad` enabled (small system for fast finite-difference gradient checking). +3. Compares values extracted from each run's `chkgrad.log`: + - Objective function value (`Function Value = ...`) + - Gradient check summary (`Maximum difference = ... Error= ...`) +4. Fails if differences exceed tolerances. + +You can override tolerances with environment variables: + +```bash +OBJ_TOL=1e-8 GRAD_TOL=1e-3 ./testing/check_numerics_profiles.sh +``` + ## References Please always cite one of the following references in publications for which Packmol was useful: @@ -93,4 +188,3 @@ JM Martinez, L Martinez, Packing optimization for the automated generation of co (https://doi.org/10.1002/jcc.10216) - diff --git a/app/packmol.f90 b/app/packmol.f90 index 2708318..df2ca32 100644 --- a/app/packmol.f90 +++ b/app/packmol.f90 @@ -593,9 +593,14 @@ program packmol keyword(iline,1).eq.'below') then nratom(icart) = nratom(icart) + 1 iratcount = iratcount + 1 + iirest = -1 do irest = 1, nrest if(irestline(irest).eq.iline) iirest = irest end do + if(iirest.eq.-1) then + write(*,*) ' ERROR: Could not map atom-level restriction to definition at line ', iline + stop exit_code_input_error + end if iratom(icart,iratcount) = iirest end if end if @@ -608,9 +613,14 @@ program packmol keyword(iline,1).eq.'below') then nratom(icart) = nratom(icart) + 1 iratcount = iratcount + 1 + iirest = -1 do irest = 1, nrest if(irestline(irest).eq.iline) iirest = irest end do + if(iirest.eq.-1) then + write(*,*) ' ERROR: Could not map restriction to definition at line ', iline + stop exit_code_input_error + end if iratom(icart,iratcount) = iirest end if end do @@ -710,6 +720,9 @@ program packmol call initial(n,x) + ! Pre-compute fixed_short_marker once (fixedatom and use_short_radius are static) + call init_fixed_short_marker() + ! Computing the energy at the initial point radscale = 1.d0 @@ -978,4 +991,3 @@ program packmol end if end program packmol - diff --git a/src/cell_indexing.f90 b/src/cell_indexing.f90 index fc8c517..7595b5a 100644 --- a/src/cell_indexing.f90 +++ b/src/cell_indexing.f90 @@ -1,7 +1,30 @@ module cell_indexing -implicit none -public :: setcell, icell_to_cell, index_cell +implicit none +public :: setcell, icell_to_cell, index_cell, n_forward_offsets, forward_offsets + +! Ordered list of relative cell offsets used for forward-only traversal. +! 1: (0, 0, 0) -> self cell (handled with latomnext to avoid double-counting) +! 2-4: forward face neighbors +! 5-10: forward edge neighbors +! 11-14: forward vertex neighbors +integer, parameter :: n_forward_offsets = 14 +integer, parameter :: forward_offsets(3,n_forward_offsets) = reshape([ & + 0, 0, 0, & + 1, 0, 0, & + 0, 1, 0, & + 0, 0, 1, & + 1,-1, 0, & + 1, 0,-1, & + 0, 1,-1, & + 0, 1, 1, & + 1, 1, 0, & + 1, 0, 1, & + 1,-1,-1, & + 1,-1, 1, & + 1, 1,-1, & + 1, 1, 1 & +], [3,n_forward_offsets]) contains ! diff --git a/src/checkpoint.f90 b/src/checkpoint.f90 index f43dff8..bdd8f47 100644 --- a/src/checkpoint.f90 +++ b/src/checkpoint.f90 @@ -12,7 +12,7 @@ subroutine checkpoint(n,x) use sizes - use compute_data + use compute_data, only : comptype, init1 use input use usegencan use ahestetic diff --git a/src/compute_data.f90 b/src/compute_data.f90 index 3bce7c0..a1c67e3 100644 --- a/src/compute_data.f90 +++ b/src/compute_data.f90 @@ -18,6 +18,7 @@ module compute_data integer, allocatable :: ityperest(:) ! (maxrest) integer, allocatable :: ibmol(:) ! (ntotat) integer, allocatable :: ibtype(:) ! (ntotat) + integer, allocatable :: fixed_short_marker(:) ! 0: movable/no-short, 1: movable/short, 2: fixed/no-short, 3: fixed/short double precision :: scale, scale2 double precision :: fdist, frest @@ -59,4 +60,12 @@ module compute_data integer, allocatable :: lcellnext(:) ! (ncells(1)*ncells(2)*ncells(3)) logical, allocatable :: empty_cell(:,:,:) ! (ncells(1),ncells(2),ncells(3)) +contains + + subroutine init_fixed_short_marker() + implicit none + if ( .not. allocated(fixed_short_marker) ) return + fixed_short_marker(:) = merge(2, 0, fixedatom(:)) + merge(1, 0, use_short_radius(:)) + end subroutine init_fixed_short_marker + end module compute_data diff --git a/src/computef.f90 b/src/computef.f90 index b5927ce..1eccdd2 100644 --- a/src/computef.f90 +++ b/src/computef.f90 @@ -9,14 +9,18 @@ subroutine computef(n,x,f) use sizes - use cell_indexing, only: index_cell, icell_to_cell, setcell - use compute_data + use cell_indexing, only: index_cell, icell_to_cell, setcell, n_forward_offsets, forward_offsets + use compute_data, only : ntotmol, ntype, ncells, nmols, natoms, idfirst, ibmol, ibtype, fdist, frest, & + xcart, coor, radius, frest_atom, comptype, init1, move, & + latomnext, latomfirst, lcellfirst, lcellnext, empty_cell use pbc implicit none integer :: n, i, j, k, icell integer :: ilugan, ilubar, icart, itype, imol, iatom, idatom - integer :: cell(3) + integer :: cell(3), neigh_cell(3) + integer :: ioffset + integer :: neigh_first(n_forward_offsets) double precision :: v1(3), v2(3), v3(3) double precision :: x(n) @@ -32,7 +36,9 @@ subroutine computef(n,x,f) ! Reset cells - if(.not.init1) call resetcells() + if(.not.init1) then + call resetcells() + end if ! Transform baricenter and angles into cartesian coordinates ! Computes cartesian coordinates from vector x and coor @@ -120,30 +126,28 @@ subroutine computef(n,x,f) j = cell(2) k = cell(3) - icart = latomfirst(i,j,k) + ! Load forward neighbor heads (self-cell + 13 forward neighbors). + do ioffset = 1, n_forward_offsets + neigh_cell(1) = cell_ind(i + forward_offsets(1,ioffset), ncells(1)) + neigh_cell(2) = cell_ind(j + forward_offsets(2,ioffset), ncells(2)) + neigh_cell(3) = cell_ind(k + forward_offsets(3,ioffset), ncells(3)) + neigh_first(ioffset) = latomfirst(neigh_cell(1),neigh_cell(2),neigh_cell(3)) + end do + + icart = neigh_first(1) do while( icart > 0 ) - if(comptype(ibtype(icart))) then - ! Interactions inside cell - f = f + fparc(icart,latomnext(icart)) - ! Interactions of cells that share faces (6 faces - 3 forward) - f = f + fparc(icart,latomfirst(cell_ind(i+1, ncells(1)),j,k)) ! 4 - (1, 0, 0) - f = f + fparc(icart,latomfirst(i,cell_ind(j+1, ncells(2)),k)) ! 5 - (0, 1, 0) - f = f + fparc(icart,latomfirst(i,j,cell_ind(k+1, ncells(3)))) ! 6 - (0, 0, 1) - ! Interactions of cells that share axes (12 edges - 6 forward) - f = f + fparc(icart,latomfirst(cell_ind(i+1, ncells(1)),cell_ind(j-1, ncells(2)),k)) ! 4 - (1, -1, 0) - f = f + fparc(icart,latomfirst(cell_ind(i+1, ncells(1)),j,cell_ind(k-1, ncells(3)))) ! 5 - (1, 0, -1) - f = f + fparc(icart,latomfirst(i,cell_ind(j+1, ncells(2)),cell_ind(k-1, ncells(3)))) ! 6 - (0, 1, -1) - f = f + fparc(icart,latomfirst(i,cell_ind(j+1, ncells(2)),cell_ind(k+1, ncells(3)))) ! 9 - (0, 1, 1) - f = f + fparc(icart,latomfirst(cell_ind(i+1, ncells(1)),cell_ind(j+1, ncells(2)),k)) ! 10 - (1, 1, 0) - f = f + fparc(icart,latomfirst(cell_ind(i+1, ncells(1)),j,cell_ind(k+1, ncells(3)))) ! 11 - (1, 0, 1) - ! Interactions of cells that share vertices (8 vertices, 8 forward) - f = f + fparc(icart,latomfirst(cell_ind(i+1, ncells(1)),cell_ind(j-1, ncells(2)),cell_ind(k-1, ncells(3)))) ! 5 - (1, -1, -1) - f = f + fparc(icart,latomfirst(cell_ind(i+1, ncells(1)),cell_ind(j-1, ncells(2)),cell_ind(k+1, ncells(3)))) ! 6 - (1, -1, 1) - f = f + fparc(icart,latomfirst(cell_ind(i+1, ncells(1)),cell_ind(j+1, ncells(2)),cell_ind(k-1, ncells(3)))) ! 7 - (1, 1, -1) - f = f + fparc(icart,latomfirst(cell_ind(i+1, ncells(1)),cell_ind(j+1, ncells(2)),cell_ind(k+1, ncells(3)))) ! 8 - (1, 1, 1) + if ( .not. comptype(ibtype(icart)) ) then + icart = latomnext(icart) + cycle end if + ! Offset #1 is self-cell and uses latomnext(icart) to keep forward-only pairs. + f = f + fparc(icart,latomnext(icart)) + do ioffset = 2, n_forward_offsets + f = f + fparc(icart,neigh_first(ioffset)) + end do + icart = latomnext(icart) end do @@ -151,4 +155,5 @@ subroutine computef(n,x,f) end do return + end subroutine computef diff --git a/src/computeg.f90 b/src/computeg.f90 index ddbddf9..1eedbf4 100644 --- a/src/computeg.f90 +++ b/src/computeg.f90 @@ -9,17 +9,20 @@ subroutine computeg(n,x,g) use sizes - use cell_indexing, only: index_cell, icell_to_cell, setcell - use compute_data + use cell_indexing, only: index_cell, icell_to_cell, setcell, n_forward_offsets, forward_offsets + use compute_data, only : ntotmol, ntype, ncells, nmols, natoms, idfirst, nratom, iratom, ibmol, ibtype, & + xcart, coor, radius, gxcar, comptype, init1, latomnext, & + latomfirst, lcellfirst, lcellnext, empty_cell use pbc implicit none integer :: n integer :: idatom, iatom, irest integer :: i, j, k, ilubar, ilugan, icart, itype, imol - integer :: icell, cell(3) + integer :: icell, cell(3), neigh_cell(3) integer :: k1, k2 - integer :: iratcount + integer :: iratcount, ioffset + integer :: neigh_first(n_forward_offsets) double precision :: x(n), g(n) double precision :: dv1beta(3), dv1gama(3), dv1teta(3),& @@ -34,7 +37,9 @@ subroutine computeg(n,x,g) ! Reset cells - if(.not.init1) call resetcells() + if(.not.init1) then + call resetcells() + end if ! Transform baricenter and angles into cartesian coordinates @@ -110,30 +115,28 @@ subroutine computeg(n,x,g) j = cell(2) k = cell(3) - icart = latomfirst(i,j,k) + ! Load forward neighbor heads (self-cell + 13 forward neighbors). + do ioffset = 1, n_forward_offsets + neigh_cell(1) = cell_ind(i + forward_offsets(1,ioffset), ncells(1)) + neigh_cell(2) = cell_ind(j + forward_offsets(2,ioffset), ncells(2)) + neigh_cell(3) = cell_ind(k + forward_offsets(3,ioffset), ncells(3)) + neigh_first(ioffset) = latomfirst(neigh_cell(1),neigh_cell(2),neigh_cell(3)) + end do + + icart = neigh_first(1) do while( icart > 0 ) - if(comptype(ibtype(icart))) then - ! Interactions inside cell - call gparc(icart,latomnext(icart)) - ! Interactions of cells that share faces (6 faces - 3 forward) - call gparc(icart,latomfirst(cell_ind(i+1, ncells(1)),j,k)) ! 4 - (1, 0, 0) - call gparc(icart,latomfirst(i,cell_ind(j+1, ncells(2)),k)) ! 5 - (0, 1, 0) - call gparc(icart,latomfirst(i,j,cell_ind(k+1, ncells(3)))) ! 6 - (0, 0, 1) - ! Interactions of cells that share axes (12 edges - 6 forward) - call gparc(icart,latomfirst(i,cell_ind(j+1, ncells(2)),cell_ind(k+1, ncells(3)))) ! 5 - (0, 1, 1) - call gparc(icart,latomfirst(i,cell_ind(j+1, ncells(2)),cell_ind(k-1, ncells(3)))) ! 6 - (0, 1, -1) - call gparc(icart,latomfirst(cell_ind(i+1, ncells(1)),cell_ind(j+1, ncells(2)),k)) ! 9 - (1, 1, 0) - call gparc(icart,latomfirst(cell_ind(i+1, ncells(1)),j,cell_ind(k+1, ncells(3)))) ! 10 - (1, 0, 1) - call gparc(icart,latomfirst(cell_ind(i+1, ncells(1)),cell_ind(j-1, ncells(2)),k)) ! 11 - (1, -1, 0) - call gparc(icart,latomfirst(cell_ind(i+1, ncells(1)),j,cell_ind(k-1, ncells(3)))) ! 12 - (1, 0, -1) - ! Interactions of cells that share vertices (8 vertices, 4 forward) - call gparc(icart,latomfirst(cell_ind(i+1, ncells(1)),cell_ind(j+1, ncells(2)),cell_ind(k+1, ncells(3)))) ! 1 - (1, 1, 1) - call gparc(icart,latomfirst(cell_ind(i+1, ncells(1)),cell_ind(j+1, ncells(2)),cell_ind(k-1, ncells(3)))) ! 2 - (1, 1, -1) - call gparc(icart,latomfirst(cell_ind(i+1, ncells(1)),cell_ind(j-1, ncells(2)),cell_ind(k+1, ncells(3)))) ! 3 - (1, -1, 1) - call gparc(icart,latomfirst(cell_ind(i+1, ncells(1)),cell_ind(j-1, ncells(2)),cell_ind(k-1, ncells(3)))) ! 4 - (1, -1, -1) + if ( .not. comptype(ibtype(icart)) ) then + icart = latomnext(icart) + cycle end if + ! Offset #1 is self-cell and uses latomnext(icart) to keep forward-only pairs. + call gparc(icart,latomnext(icart)) + do ioffset = 2, n_forward_offsets + call gparc(icart,neigh_first(ioffset)) + end do + icart = latomnext(icart) end do @@ -241,5 +244,5 @@ subroutine computeg(n,x,g) end do return -end subroutine computeg +end subroutine computeg diff --git a/src/fparc.f90 b/src/fparc.f90 index 13429f1..addff94 100644 --- a/src/fparc.f90 +++ b/src/fparc.f90 @@ -10,8 +10,9 @@ double precision function fparc(icart,firstjcart) use sizes - use compute_data - use pbc, only : delta_vector, pbc_length + use compute_data, only : comptype, fdist, fdist_atom, move, radius, radius_ini, fscale, short_radius, & + short_radius_scale, fixed_short_marker, latomnext, xcart, ibmol, ibtype + use pbc, only : pbc_length implicit none ! SCALAR ARGUMENTS @@ -20,59 +21,177 @@ double precision function fparc(icart,firstjcart) ! LOCAL SCALARS integer :: jcart double precision :: datom, tol, short_tol, short_tol_penalty, short_tol_scale - double precision :: vdiff(3) + double precision :: sumr, short_sumr, delta, tol_sq, short_tol_sq + double precision :: dx, dy, dz + double precision :: xi, yi, zi, ri, rii, fsi, short_ri, short_si + double precision :: inv_pbcx, inv_pbcy, inv_pbcz + integer :: marker_i, itype_i, imol_i + logical :: use_short_i, fixed_i + + xi = xcart(icart,1) + yi = xcart(icart,2) + zi = xcart(icart,3) + ri = radius(icart) + rii = radius_ini(icart) + fsi = fscale(icart) + short_ri = short_radius(icart) + short_si = short_radius_scale(icart) + marker_i = fixed_short_marker(icart) + use_short_i = iand(marker_i,1) /= 0 + fixed_i = marker_i >= 2 + itype_i = ibtype(icart) + imol_i = ibmol(icart) + inv_pbcx = 1.d0 / pbc_length(1) + inv_pbcy = 1.d0 / pbc_length(2) + inv_pbcz = 1.d0 / pbc_length(3) fparc = 0.0d0 - jcart = firstjcart - do while ( jcart > 0 ) - ! - ! Cycle if this type is not to be computed - ! - if ( .not. comptype(ibtype(jcart))) then - jcart = latomnext(jcart) - cycle - end if - ! - ! Cycle if the atoms are from the same molecule - ! - if ( ibmol(icart) == ibmol(jcart) .and. & - ibtype(icart) == ibtype(jcart) ) then - jcart = latomnext(jcart) - cycle - end if - ! - ! Cycle if both atoms are from fixed molecules - ! - if ( fixedatom(icart) .and. fixedatom(jcart) ) then + if ( .not. fixed_i .and. .not. use_short_i ) then + ! Fast path: movable atom i and no short-radius term from i. + jcart = firstjcart + do while ( jcart > 0 ) + if ( .not. comptype(ibtype(jcart))) then + jcart = latomnext(jcart) + cycle + end if + if ( imol_i == ibmol(jcart) .and. itype_i == ibtype(jcart) ) then + jcart = latomnext(jcart) + cycle + end if + + dx = xi - xcart(jcart,1) + dy = yi - xcart(jcart,2) + dz = zi - xcart(jcart,3) + dx = dx - pbc_length(1) * dnint(dx * inv_pbcx) + dy = dy - pbc_length(2) * dnint(dy * inv_pbcy) + dz = dz - pbc_length(3) * dnint(dz * inv_pbcz) + datom = dx*dx + dy*dy + dz*dz + + sumr = ri + radius(jcart) + tol = sumr*sumr + if ( datom < tol ) then + delta = datom - tol + fparc = fparc + fsi*fscale(jcart)*delta*delta + if ( iand(fixed_short_marker(jcart),1) /= 0 ) then + short_sumr = short_ri + short_radius(jcart) + short_tol = short_sumr*short_sumr + if ( datom < short_tol ) then + short_tol_penalty = datom-short_tol + short_tol_scale = dsqrt(short_si*short_radius_scale(jcart)) + tol_sq = tol*tol + short_tol_sq = short_tol*short_tol + short_tol_scale = short_tol_scale*(tol_sq/short_tol_sq) + fparc = fparc + fsi*fscale(jcart)*short_tol_scale*short_tol_penalty*short_tol_penalty + end if + end if + end if + sumr = rii + radius_ini(jcart) + tol = sumr*sumr + fdist = dmax1(tol-datom,fdist) + if ( move ) then + fdist_atom(icart) = dmax1(tol-datom,fdist_atom(icart)) + fdist_atom(jcart) = dmax1(tol-datom,fdist_atom(jcart)) + end if jcart = latomnext(jcart) - cycle - end if - ! - ! Otherwise, compute distance and evaluate function for this pair - ! - vdiff = delta_vector(xcart(icart,:), xcart(jcart,:), pbc_length) - datom = ( vdiff(1) )**2 + ( vdiff(2) )**2 + ( vdiff(3) )**2 - tol = (radius(icart)+radius(jcart))**2 - if ( datom < tol ) then - fparc = fparc + fscale(icart)*fscale(jcart)*(datom-tol)**2 - if ( use_short_radius(icart) .or. use_short_radius(jcart) ) then - short_tol = (short_radius(icart)+short_radius(jcart))**2 + end do + else if ( .not. fixed_i ) then + ! Short-radius-enabled path for movable atoms. + jcart = firstjcart + do while ( jcart > 0 ) + if ( .not. comptype(ibtype(jcart))) then + jcart = latomnext(jcart) + cycle + end if + if ( imol_i == ibmol(jcart) .and. itype_i == ibtype(jcart) ) then + jcart = latomnext(jcart) + cycle + end if + + dx = xi - xcart(jcart,1) + dy = yi - xcart(jcart,2) + dz = zi - xcart(jcart,3) + dx = dx - pbc_length(1) * dnint(dx * inv_pbcx) + dy = dy - pbc_length(2) * dnint(dy * inv_pbcy) + dz = dz - pbc_length(3) * dnint(dz * inv_pbcz) + datom = dx*dx + dy*dy + dz*dz + + sumr = ri + radius(jcart) + tol = sumr*sumr + if ( datom < tol ) then + delta = datom - tol + fparc = fparc + fsi*fscale(jcart)*delta*delta + short_sumr = short_ri + short_radius(jcart) + short_tol = short_sumr*short_sumr if ( datom < short_tol ) then short_tol_penalty = datom-short_tol - short_tol_scale = dsqrt(short_radius_scale(icart)*short_radius_scale(jcart)) - short_tol_scale = short_tol_scale*(tol**2/short_tol**2) - fparc = fparc + fscale(icart)*fscale(jcart)*short_tol_scale*short_tol_penalty**2 + short_tol_scale = dsqrt(short_si*short_radius_scale(jcart)) + tol_sq = tol*tol + short_tol_sq = short_tol*short_tol + short_tol_scale = short_tol_scale*(tol_sq/short_tol_sq) + fparc = fparc + fsi*fscale(jcart)*short_tol_scale*short_tol_penalty*short_tol_penalty end if end if - end if - tol = (radius_ini(icart)+radius_ini(jcart))**2 - fdist = dmax1(tol-datom,fdist) - if ( move ) then - fdist_atom(icart) = dmax1(tol-datom,fdist_atom(icart)) - fdist_atom(jcart) = dmax1(tol-datom,fdist_atom(jcart)) - end if - jcart = latomnext(jcart) - end do + sumr = rii + radius_ini(jcart) + tol = sumr*sumr + fdist = dmax1(tol-datom,fdist) + if ( move ) then + fdist_atom(icart) = dmax1(tol-datom,fdist_atom(icart)) + fdist_atom(jcart) = dmax1(tol-datom,fdist_atom(jcart)) + end if + jcart = latomnext(jcart) + end do + else + ! Fixed-atom path (keeps explicit fixed/fixed exclusion). + jcart = firstjcart + do while ( jcart > 0 ) + if ( .not. comptype(ibtype(jcart))) then + jcart = latomnext(jcart) + cycle + end if + if ( imol_i == ibmol(jcart) .and. itype_i == ibtype(jcart) ) then + jcart = latomnext(jcart) + cycle + end if + if ( fixed_short_marker(jcart) >= 2 ) then + jcart = latomnext(jcart) + cycle + end if -end function fparc + dx = xi - xcart(jcart,1) + dy = yi - xcart(jcart,2) + dz = zi - xcart(jcart,3) + dx = dx - pbc_length(1) * dnint(dx * inv_pbcx) + dy = dy - pbc_length(2) * dnint(dy * inv_pbcy) + dz = dz - pbc_length(3) * dnint(dz * inv_pbcz) + datom = dx*dx + dy*dy + dz*dz + sumr = ri + radius(jcart) + tol = sumr*sumr + if ( datom < tol ) then + delta = datom - tol + fparc = fparc + fsi*fscale(jcart)*delta*delta + if ( use_short_i .or. iand(fixed_short_marker(jcart),1) /= 0 ) then + short_sumr = short_ri + short_radius(jcart) + short_tol = short_sumr*short_sumr + if ( datom < short_tol ) then + short_tol_penalty = datom-short_tol + short_tol_scale = dsqrt(short_si*short_radius_scale(jcart)) + tol_sq = tol*tol + short_tol_sq = short_tol*short_tol + short_tol_scale = short_tol_scale*(tol_sq/short_tol_sq) + fparc = fparc + fsi*fscale(jcart)*short_tol_scale*short_tol_penalty*short_tol_penalty + end if + end if + end if + sumr = rii + radius_ini(jcart) + tol = sumr*sumr + fdist = dmax1(tol-datom,fdist) + if ( move ) then + fdist_atom(icart) = dmax1(tol-datom,fdist_atom(icart)) + fdist_atom(jcart) = dmax1(tol-datom,fdist_atom(jcart)) + end if + jcart = latomnext(jcart) + end do + end if + +end function fparc diff --git a/src/gencan.f b/src/gencan.f index 389a284..1c514f3 100644 --- a/src/gencan.f +++ b/src/gencan.f @@ -49,9 +49,16 @@ subroutine evalal(n,x,m,lambda,rho,f,flag) double precision lambda(m),rho(m),x(n) C LOCAL SCALARS + double precision ignore flag = 0 +C Keep ALGENCAN-compatible dummy arguments intentionally referenced. + ignore = 0.0d0 + if ( m .gt. 0 ) then + ignore = ignore + lambda(1) - lambda(1) + rho(1) - rho(1) + end if + call computef(n,x,f) end @@ -107,9 +114,16 @@ subroutine evalnal(n,x,m,lambda,rho,g,flag) double precision g(n),lambda(m),rho(m),x(n) C LOCAL SCALARS + double precision ignore flag = 0 +C Keep ALGENCAN-compatible dummy arguments intentionally referenced. + ignore = 0.0d0 + if ( m .gt. 0 ) then + ignore = ignore + lambda(1) - lambda(1) + rho(1) - rho(1) + end if + call computeg(n,x,g) end @@ -117,12 +131,7 @@ subroutine evalnal(n,x,m,lambda,rho,g,flag) C ***************************************************************** C ***************************************************************** -c Modified by L. Martinez (there was an error on the number of -c parameters when calling this subroutine). This subroutine does -c nothing. -c subroutine evalhd(nind,ind,n,x,m,lambda,rho,d,hd,flag) - - subroutine evalhd(n) + subroutine evalhd(nind,ind,n,x,m,lambda,rho,d,hd,flag) C This subroutine computes the product of the Hessian matrix times C the input vector argument d. If GENCAN is being used stand-alone @@ -184,14 +193,28 @@ subroutine evalhd(n) implicit none C SCALAR ARGUMENTS -c integer flag,m,n,nind - integer n + integer flag,m,n,nind C ARRAY ARGUMENTS -c integer ind(nind) -c double precision d(n),hd(n),lambda(m),rho(m),x(n) + integer ind(nind) + double precision d(n),hd(n),lambda(m),rho(m),x(n) + +C LOCAL SCALARS + double precision ignore -c flag = - 1 +C Keep the full evalhd API for compatibility with GENCAN/ALGENCAN. +C This default stub intentionally does not compute Hessian products. + flag = 0 + ignore = 0.0d0 + if ( nind .gt. 0 ) then + ignore = ignore + ind(1) - ind(1) + end if + if ( n .gt. 0 ) then + ignore = ignore + d(1) - d(1) + x(1) - x(1) + hd(1) - hd(1) + end if + if ( m .gt. 0 ) then + ignore = ignore + lambda(1) - lambda(1) + rho(1) - rho(1) + end if end @@ -1719,6 +1742,11 @@ subroutine gencan(n,x,l,u,m,lambda,rho,epsgpen,epsgpsn,maxitnfp, bestprog = 0.0d0 itnfp = 0 +C Initialize BB-related scalars. They are updated at the end of each +C iteration; zero defaults force first-iteration fallback formulas. + sts = 0.0d0 + sty = 0.0d0 + C for testing progress in the projected gradient norm. do i = 0,maxitngp - 1 lastgpns(i) = infabs @@ -2125,7 +2153,9 @@ subroutine gencan(n,x,l,u,m,lambda,rho,epsgpen,epsgpsn,maxitnfp, delta = udelta0 end if else - delta = max( delmin, 10.0d0 * sqrt( sts ) ) +C sts= should be non-negative; guard protects against +C accidental first-use reads and round-off negatives. + delta = max( delmin, 10.0d0 * sqrt( max( 0.0d0, sts ) ) ) end if C Shrink the point, its gradient and the bounds @@ -3244,6 +3274,18 @@ subroutine cg(nind,ind,n,x,m,lambda,rho,g,delta,l,u,eps,epsnqmp, snorm2 = 0.0d0 rnorm2 = gnorm2 +C Safe defaults for quantities that are only meaningful after the +C first CG step update. These values preserve first-iteration logic. + alpha = 0.0d0 + dnorm2 = 0.0d0 + dtr = 0.0d0 + dtw = 0.0d0 + rnorm2prev = rnorm2 + rbdposaind = 1 + rbdposatype = 1 + rbdnegaind = 1 + rbdnegatype = 1 + C ================================================================== C Print initial information C ================================================================== @@ -3329,7 +3371,13 @@ subroutine cg(nind,ind,n,x,m,lambda,rho,g,delta,l,u,eps,epsnqmp, else - beta = rnorm2 / rnorm2prev +C Polak-Ribiere-like coefficient requires previous residual norm. +C If it is unavailable/non-positive, restart with beta=0. + if ( rnorm2prev .gt. 0.0d0 ) then + beta = rnorm2 / rnorm2prev + else + beta = 0.0d0 + end if do i = 1,nind d(i) = - r(i) + beta * d(i) @@ -3625,6 +3673,8 @@ subroutine cg(nind,ind,n,x,m,lambda,rho,g,delta,l,u,eps,epsnqmp, if ( alpha .eq. amax2 .or. alpha .eq. amax2n ) then +C Boundary metadata is set when a finite box step is identified. +C Keep initialized fallback values for degenerate directions. if ( alpha .eq. amax2 ) then rbdind = rbdposaind rbdtype = rbdposatype @@ -5214,7 +5264,7 @@ subroutine calchd(nind,ind,x,d,g,n,xc,m,lambda,rho,hd,xtmp,sterel, C Compute the Hessian times vector d product calling the user C supplied subroutine evalhd - call evalhd(n) + call evalhd(nind,ind,n,x,m,lambda,rho,d,hd,inform) C Shrink x, d and hd to the reduced space @@ -5989,4 +6039,4 @@ DOUBLE PRECISION FUNCTION HSLDNRM2(N,DX,INCX) C near to the solution and less strict when we ar far from the C solution. c -C 8) Many things in the output were changed. \ No newline at end of file +C 8) Many things in the output were changed. diff --git a/src/gparc.f90 b/src/gparc.f90 index 9a68a76..54948b2 100644 --- a/src/gparc.f90 +++ b/src/gparc.f90 @@ -9,8 +9,9 @@ subroutine gparc(icart,firstjcart) use sizes - use compute_data - use pbc, only : delta_vector, pbc_length + use compute_data, only : comptype, radius, fscale, short_radius, short_radius_scale, fixed_short_marker, & + latomnext, gxcar, xcart, ibmol, ibtype + use pbc, only : pbc_length implicit none ! SCALAR ARGUMENTS @@ -19,69 +20,203 @@ subroutine gparc(icart,firstjcart) ! LOCAL SCALARS integer :: jcart double precision :: datom, dtemp, xdiff, tol, short_tol, short_tol_scale - double precision :: vdiff(3) + double precision :: sumr, short_sumr, tol_sq, short_tol_sq + double precision :: dx, dy, dz + double precision :: xi, yi, zi, ri, fsi, short_ri, short_si + double precision :: inv_pbcx, inv_pbcy, inv_pbcz + integer :: marker_i, itype_i, imol_i + logical :: use_short_i, fixed_i + + xi = xcart(icart,1) + yi = xcart(icart,2) + zi = xcart(icart,3) + ri = radius(icart) + fsi = fscale(icart) + short_ri = short_radius(icart) + short_si = short_radius_scale(icart) + marker_i = fixed_short_marker(icart) + use_short_i = iand(marker_i,1) /= 0 + fixed_i = marker_i >= 2 + itype_i = ibtype(icart) + imol_i = ibmol(icart) + inv_pbcx = 1.d0 / pbc_length(1) + inv_pbcy = 1.d0 / pbc_length(2) + inv_pbcz = 1.d0 / pbc_length(3) + + if ( .not. fixed_i .and. .not. use_short_i ) then + ! Fast path: movable atom i and no short-radius term from i. + jcart = firstjcart + do while ( jcart > 0 ) + if ( .not. comptype(ibtype(jcart))) then + jcart = latomnext(jcart) + cycle + end if + if ( imol_i == ibmol(jcart) .and. itype_i == ibtype(jcart) ) then + jcart = latomnext(jcart) + cycle + end if + + sumr = ri + radius(jcart) + tol = sumr*sumr + dx = xi - xcart(jcart,1) + dy = yi - xcart(jcart,2) + dz = zi - xcart(jcart,3) + dx = dx - pbc_length(1) * dnint(dx * inv_pbcx) + dy = dy - pbc_length(2) * dnint(dy * inv_pbcy) + dz = dz - pbc_length(3) * dnint(dz * inv_pbcz) + datom = dx*dx + dy*dy + dz*dz + if( datom < tol ) then + dtemp = fsi*fscale(jcart) * 4.d0 * (datom - tol) + xdiff = dtemp*dx + gxcar(icart,1)= gxcar(icart,1) + xdiff + gxcar(jcart,1)= gxcar(jcart,1) - xdiff + xdiff = dtemp*dy + gxcar(icart,2)= gxcar(icart,2) + xdiff + gxcar(jcart,2)= gxcar(jcart,2) - xdiff + xdiff = dtemp*dz + gxcar(icart,3)= gxcar(icart,3) + xdiff + gxcar(jcart,3)= gxcar(jcart,3) - xdiff + if ( iand(fixed_short_marker(jcart),1) /= 0 ) then + short_sumr = short_ri + short_radius(jcart) + short_tol = short_sumr*short_sumr + if ( datom < short_tol ) then + short_tol_scale = dsqrt(short_si*short_radius_scale(jcart)) + tol_sq = tol*tol + short_tol_sq = short_tol*short_tol + short_tol_scale = short_tol_scale*( tol_sq / short_tol_sq ) + dtemp = fsi*fscale(jcart) * short_tol_scale * 4.d0 * (datom - short_tol) + xdiff = dtemp*dx + gxcar(icart,1)= gxcar(icart,1) + xdiff + gxcar(jcart,1)= gxcar(jcart,1) - xdiff + xdiff = dtemp*dy + gxcar(icart,2)= gxcar(icart,2) + xdiff + gxcar(jcart,2)= gxcar(jcart,2) - xdiff + xdiff = dtemp*dz + gxcar(icart,3)= gxcar(icart,3) + xdiff + gxcar(jcart,3)= gxcar(jcart,3) - xdiff + end if + end if + end if - jcart = firstjcart - do while ( jcart > 0 ) - ! - ! Cycle if this type is not to be computed - ! - if ( .not. comptype(ibtype(jcart))) then - jcart = latomnext(jcart) - cycle - end if - ! - ! Cycle if the atoms are from the same molecule - ! - if ( ibmol(icart) == ibmol(jcart) .and. & - ibtype(icart) == ibtype(jcart) ) then - jcart = latomnext(jcart) - cycle - end if - ! - ! Cycle if both atoms are from fixed molecules - ! - if ( fixedatom(icart) .and. fixedatom(jcart) ) then jcart = latomnext(jcart) - cycle - end if - ! - ! Otherwise, compute distance and evaluate function for this pair - ! - tol = (radius(icart)+radius(jcart))**2 - vdiff = delta_vector(xcart(icart,:), xcart(jcart,:), pbc_length) - datom = (vdiff(1))**2 + (vdiff(2))**2 + (vdiff(3))**2 - if( datom < tol ) then - dtemp = fscale(icart)*fscale(jcart) * 4.d0 * (datom - tol) - xdiff = dtemp*vdiff(1) - gxcar(icart,1)= gxcar(icart,1) + xdiff - gxcar(jcart,1)= gxcar(jcart,1) - xdiff - xdiff = dtemp*vdiff(2) - gxcar(icart,2)= gxcar(icart,2) + xdiff - gxcar(jcart,2)= gxcar(jcart,2) - xdiff - xdiff = dtemp*vdiff(3) - gxcar(icart,3)= gxcar(icart,3) + xdiff - gxcar(jcart,3)= gxcar(jcart,3) - xdiff - if ( use_short_radius(icart) .or. use_short_radius(jcart) ) then - short_tol = ( short_radius(icart) + short_radius(jcart) )**2 + end do + else if ( .not. fixed_i ) then + ! Short-radius-enabled path for movable atoms. + jcart = firstjcart + do while ( jcart > 0 ) + if ( .not. comptype(ibtype(jcart))) then + jcart = latomnext(jcart) + cycle + end if + if ( imol_i == ibmol(jcart) .and. itype_i == ibtype(jcart) ) then + jcart = latomnext(jcart) + cycle + end if + + sumr = ri + radius(jcart) + tol = sumr*sumr + dx = xi - xcart(jcart,1) + dy = yi - xcart(jcart,2) + dz = zi - xcart(jcart,3) + dx = dx - pbc_length(1) * dnint(dx * inv_pbcx) + dy = dy - pbc_length(2) * dnint(dy * inv_pbcy) + dz = dz - pbc_length(3) * dnint(dz * inv_pbcz) + datom = dx*dx + dy*dy + dz*dz + if( datom < tol ) then + dtemp = fsi*fscale(jcart) * 4.d0 * (datom - tol) + xdiff = dtemp*dx + gxcar(icart,1)= gxcar(icart,1) + xdiff + gxcar(jcart,1)= gxcar(jcart,1) - xdiff + xdiff = dtemp*dy + gxcar(icart,2)= gxcar(icart,2) + xdiff + gxcar(jcart,2)= gxcar(jcart,2) - xdiff + xdiff = dtemp*dz + gxcar(icart,3)= gxcar(icart,3) + xdiff + gxcar(jcart,3)= gxcar(jcart,3) - xdiff + short_sumr = short_ri + short_radius(jcart) + short_tol = short_sumr*short_sumr if ( datom < short_tol ) then - short_tol_scale = dsqrt(short_radius_scale(icart)*short_radius_scale(jcart)) - short_tol_scale = short_tol_scale*( tol**2 / short_tol**2 ) - dtemp = fscale(icart)*fscale(jcart) * 4.d0 * short_tol_scale*(datom - short_tol) - xdiff = dtemp*vdiff(1) + short_tol_scale = dsqrt(short_si*short_radius_scale(jcart)) + tol_sq = tol*tol + short_tol_sq = short_tol*short_tol + short_tol_scale = short_tol_scale*( tol_sq / short_tol_sq ) + dtemp = fsi*fscale(jcart) * short_tol_scale * 4.d0 * (datom - short_tol) + xdiff = dtemp*dx gxcar(icart,1)= gxcar(icart,1) + xdiff gxcar(jcart,1)= gxcar(jcart,1) - xdiff - xdiff = dtemp*vdiff(2) + xdiff = dtemp*dy gxcar(icart,2)= gxcar(icart,2) + xdiff gxcar(jcart,2)= gxcar(jcart,2) - xdiff - xdiff = dtemp*vdiff(3) + xdiff = dtemp*dz gxcar(icart,3)= gxcar(icart,3) + xdiff gxcar(jcart,3)= gxcar(jcart,3) - xdiff end if end if - end if - jcart = latomnext(jcart) - end do + + jcart = latomnext(jcart) + end do + else + ! Fixed-atom path (keeps explicit fixed/fixed exclusion). + jcart = firstjcart + do while ( jcart > 0 ) + if ( .not. comptype(ibtype(jcart))) then + jcart = latomnext(jcart) + cycle + end if + if ( imol_i == ibmol(jcart) .and. itype_i == ibtype(jcart) ) then + jcart = latomnext(jcart) + cycle + end if + if ( fixed_short_marker(jcart) >= 2 ) then + jcart = latomnext(jcart) + cycle + end if + + sumr = ri + radius(jcart) + tol = sumr*sumr + dx = xi - xcart(jcart,1) + dy = yi - xcart(jcart,2) + dz = zi - xcart(jcart,3) + dx = dx - pbc_length(1) * dnint(dx * inv_pbcx) + dy = dy - pbc_length(2) * dnint(dy * inv_pbcy) + dz = dz - pbc_length(3) * dnint(dz * inv_pbcz) + datom = dx*dx + dy*dy + dz*dz + if( datom < tol ) then + dtemp = fsi*fscale(jcart) * 4.d0 * (datom - tol) + xdiff = dtemp*dx + gxcar(icart,1)= gxcar(icart,1) + xdiff + gxcar(jcart,1)= gxcar(jcart,1) - xdiff + xdiff = dtemp*dy + gxcar(icart,2)= gxcar(icart,2) + xdiff + gxcar(jcart,2)= gxcar(jcart,2) - xdiff + xdiff = dtemp*dz + gxcar(icart,3)= gxcar(icart,3) + xdiff + gxcar(jcart,3)= gxcar(jcart,3) - xdiff + if ( use_short_i .or. iand(fixed_short_marker(jcart),1) /= 0 ) then + short_sumr = short_ri + short_radius(jcart) + short_tol = short_sumr*short_sumr + if ( datom < short_tol ) then + short_tol_scale = dsqrt(short_si*short_radius_scale(jcart)) + tol_sq = tol*tol + short_tol_sq = short_tol*short_tol + short_tol_scale = short_tol_scale*( tol_sq / short_tol_sq ) + dtemp = fsi*fscale(jcart) * short_tol_scale * 4.d0 * (datom - short_tol) + xdiff = dtemp*dx + gxcar(icart,1)= gxcar(icart,1) + xdiff + gxcar(jcart,1)= gxcar(jcart,1) - xdiff + xdiff = dtemp*dy + gxcar(icart,2)= gxcar(icart,2) + xdiff + gxcar(jcart,2)= gxcar(jcart,2) - xdiff + xdiff = dtemp*dz + gxcar(icart,3)= gxcar(icart,3) + xdiff + gxcar(jcart,3)= gxcar(jcart,3) - xdiff + end if + end if + end if + + jcart = latomnext(jcart) + end do + end if + return end subroutine gparc - diff --git a/src/gwalls.f90 b/src/gwalls.f90 index 2ca127a..455330c 100644 --- a/src/gwalls.f90 +++ b/src/gwalls.f90 @@ -9,7 +9,7 @@ subroutine gwalls(icart,irest) use sizes - use compute_data + use compute_data, only : ityperest, restpars, xcart, gxcar, scale, scale2 implicit none integer :: icart, irest diff --git a/src/heuristics.f90 b/src/heuristics.f90 index d202040..a9e0f0b 100644 --- a/src/heuristics.f90 +++ b/src/heuristics.f90 @@ -9,7 +9,8 @@ subroutine movebad(n,x,fx,movebadprint) use sizes - use compute_data + use compute_data, only : ntotmol, ntype, ntotat, nmols, natoms, sizemin, sizemax, radius, radius_ini, & + fdist_atom, frest_atom, dmax, comptype, move, fmol, radiuswork use input, only : movefrac, movebadrandom, precision, maxmove use usegencan use flashsort diff --git a/src/initial.f90 b/src/initial.f90 index 2f8237e..76c8d2f 100644 --- a/src/initial.f90 +++ b/src/initial.f90 @@ -1,558 +1,563 @@ -! -! Written by Leandro Martínez, 2009-2011. -! Copyright (c) 2009-2018, Leandro Martínez, Jose Mario Martinez, -! Ernesto G. Birgin. -! -! Subroutine initial: Subroutine that reset parameters and -! builds the initial point -! - -subroutine initial(n,x) - - use exit_codes - use sizes - use cell_indexing, only: setcell - use compute_data - use input, only : randini, ntype_with_fixed, fix, chkgrad, avoidoverlap,& - discale, precision, sidemax, restart_from, input_itype,& - nloop0_type - use usegencan - use ahestetic - use pbc - implicit none - integer :: n, i, j, idatom, iatom, ilubar, ilugan, icart, itype, & - imol, ntry, cell(3), ic, jc, kc, ifatom, & - idfatom, iftype, jatom, ioerr, max_guess_try - - double precision :: x(n), beta, gamma, theta, & - fx, cell_side, rnd, xrnd(3), v1(3), v2(3), v3(3), xcm(3) - double precision :: scale_rnd_guess - double precision, parameter :: twopi = 8.d0*datan(1.d0) - double precision, allocatable :: cm_min(:,:), cm_max(:,:) - - logical :: overlap, movebadprint, hasbad - - character(len=strl) :: record - - ! We need to initialize the move logical variable - move = .false. - - ! Default status of the function evaluation - init1 = .false. - lcellfirst = 0 - - ! Initialize the comptype logical array - comptype(1:ntype_with_fixed) = .true. - - ! Penalty factors for the objective function relative to restrictions - ! Default values: scale = 1.d2, scale2 = 1.d1 - - scale = 1.d0 - scale2 = 1.d-2 - call tobar() - - ! Compute maximum internal distance within each type of molecule - - do itype = 1, ntype - dmax(itype) = 0.d0 - idatom = idfirst(itype) - 1 - do iatom = 1, natoms(itype) - 1 - do jatom = iatom + 1, natoms(itype) - dmax(itype) = dmax1 ( dmax(itype),& - (coor(idatom+iatom,1)-coor(idatom+jatom,1))**2+& - (coor(idatom+iatom,2)-coor(idatom+jatom,2))**2+& - (coor(idatom+iatom,3)-coor(idatom+jatom,3))**2 ) - end do - end do - dmax(itype) = dsqrt(dmax(itype)) - write(*,*) ' Maximum internal distance of type ',itype,': ',& - dmax(itype) - if(dmax(itype).eq.0.) dmax(itype) = 1.d0 - end do - - ! Maximum size of the system: if you system is very large (about - ! 80 nm wide), increase the sidemax parameter. - ! Otherwise, the packing can be slow and unsucesful - allocate(cm_min(ntype,3), cm_max(ntype,3)) - x(:) = 0.d0 - call restmol(1,0,n,x,fx,.true.) - sizemin(1:3) = x(1:3) - sidemax - sizemax(1:3) = x(1:3) + sidemax - write(*,*) ' All atoms must be within these coordinates: ' - write(*,*) ' x: [ ', sizemin(1),', ', sizemax(1), ' ] ' - write(*,*) ' y: [ ', sizemin(2),', ', sizemax(2), ' ] ' - write(*,*) ' z: [ ', sizemin(3),', ', sizemax(3), ' ] ' - write(*,*) ' If the system is larger than this, increase the sidemax parameter. ' - - ! Create first aleatory guess - - i = 0 - j = ntotmol*3 - do itype = 1, ntype - do imol = 1, nmols(itype) - x(i+1) = sizemin(1) + rnd()*(sizemax(1)-sizemin(1)) - x(i+2) = sizemin(2) + rnd()*(sizemax(2)-sizemin(2)) - x(i+3) = sizemin(3) + rnd()*(sizemax(3)-sizemin(3)) - if ( constrain_rot(itype,1) ) then - x(j+1) = ( rot_bound(itype,1,1) - dabs(rot_bound(itype,1,2)) ) + & - 2.d0*rnd()*dabs(rot_bound(itype,1,2)) - else - x(j+1) = twopi*rnd() - end if - if ( constrain_rot(itype,2) ) then - x(j+2) = ( rot_bound(itype,2,1) - dabs(rot_bound(itype,2,2)) ) + & - 2.d0*rnd()*dabs(rot_bound(itype,2,2)) - else - x(j+2) = twopi*rnd() - end if - if ( constrain_rot(itype,3) ) then - x(j+3) = ( rot_bound(itype,3,1) - dabs(rot_bound(itype,3,2)) ) + & - 2.d0*rnd()*dabs(rot_bound(itype,3,2)) - else - x(j+3) = twopi*rnd() - end if - i = i + 3 - j = j + 3 - end do - end do - - ! Initialize cartesian coordinate array for the first time - - ilubar = 0 - ilugan = ntotmol*3 - icart = 0 - do itype = 1, ntype - do imol = 1, nmols(itype) - xcm = x(ilubar+1:ilubar+3) - beta = x(ilugan+1) - gamma = x(ilugan+2) - theta = x(ilugan+3) - call eulerrmat(beta,gamma,theta,v1,v2,v3) - idatom = idfirst(itype) - 1 - do iatom = 1, natoms(itype) - icart = icart + 1 - idatom = idatom + 1 - call compcart(xcart(icart,1:3),xcm,coor(idatom,1:3),v1,v2,v3) - fixedatom(icart) = .false. - end do - end do - end do - if(fix) then - icart = ntotat - nfixedat - do iftype = ntype + 1, ntype_with_fixed - idfatom = idfirst(iftype) - 1 - do ifatom = 1, natoms(iftype) - idfatom = idfatom + 1 - icart = icart + 1 - xcart(icart,1) = coor(idfatom,1) - xcart(icart,2) = coor(idfatom,2) - xcart(icart,3) = coor(idfatom,3) - fixedatom(icart) = .true. - ! Check if fixed molecules are compatible with PBC given - if (using_pbc) then - do i = 1, 3 - if (xcart(icart, i) < pbc_min(i) .or. xcart(icart, i) > pbc_max(i)) then - write(*,*) "ERROR: Fixed molecule are outside the PBC box:" - write(*,*) " Atom: ", ifatom, " of molecule: ", input_itype(iftype), " - coordinate: ", i - write(*,*) " ", xcart(icart, i), " not in [", pbc_min(i), ", ", pbc_max(i), "]" - write(*,*) "(after translating/rotation the fixed molecule with the given orientation)" - stop exit_code_input_error - end if - end do - end if - end do - end do - end if - - ! Performing some steps of optimization for the restrictions only - - write(*,hash3_line) - write(*,"(' Building initial approximation ... ' )") - write(*,hash3_line) - write(*,"(' Adjusting initial point to fit the constraints ')") - write(*,dash2_line) - init1 = .true. - call swaptype(n,x,itype,0) ! Initialize swap arrays - itype = 0 - do while( itype <= ntype ) - itype = itype + 1 - if ( itype <= ntype ) then - call swaptype(n,x,itype,1) ! Set arrays for this type - else - call swaptype(n,x,itype,3) ! Restore arrays if itype = ntype + 1 - exit - end if - write(*,dash3_line) - write(*,*) ' Molecules of type: ', input_itype(itype) - write(*,*) - i = 0 - hasbad = .true. - call computef(n,x,fx) - do while( frest > precision .and. i.le. nloop0_type(itype)-1 .and. hasbad) - i = i + 1 - write(*,prog1_line) - call pgencan(n,x,fx) - call computef(n,x,fx) - if(frest > precision) then - write(*,"( a,i6,a,i6 )")' Fixing bad orientations ... ', i,' of ', nloop0_type(itype) - movebadprint = .true. - call movebad(n,x,fx,movebadprint) - end if - end do - write(*,*) - write(*,*) ' Restraint-only function value: ', fx - write(*,*) ' Maximum violation of the restraints: ', frest - - call swaptype(n,x,itype,2) ! Save current type results - - if( hasbad .and. frest > precision ) then - write(*,*) ' ERROR: Packmol was unable to put the molecules' - write(*,*) ' in the desired regions even without' - write(*,*) ' considering distance tolerances. ' - write(*,*) ' Probably there is something wrong with' - write(*,*) ' the constraints, since it seems that' - write(*,*) ' the molecules cannot satisfy them at' - write(*,*) ' at all. ' - write(*,*) ' Please check the spatial constraints and' - write(*,*) ' try again.' - if ( i .ge. nloop0_type(itype)-1 ) then - end if - write(*,*) ' >The maximum number of cycles (',nloop0_type(itype),') was achieved.' - write(*,*) ' You may try increasing it with the',' nloop0 keyword, as in: nloop0 1000 ' - stop exit_code_failed_to_converge - end if - end do - init1 = .false. - - ! Rescaling sizemin and sizemax in order to build the patch of cells - - write(*,dash3_line) - write(*,*) ' Rescaling maximum and minimum coordinates... ' - sizemin(1:3) = 1.d20 - sizemax(1:3) = -1.d20 - ! Maximum and minimum coordinates of fixed molecules - icart = ntotat - nfixedat - do itype = ntype + 1, ntype_with_fixed - do imol = 1, nmols(itype) - do iatom = 1, natoms(itype) - icart = icart + 1 - sizemin(1) = dmin1(sizemin(1),xcart(icart,1)) - sizemin(2) = dmin1(sizemin(2),xcart(icart,2)) - sizemin(3) = dmin1(sizemin(3),xcart(icart,3)) - sizemax(1) = dmax1(sizemax(1),xcart(icart,1)) - sizemax(2) = dmax1(sizemax(2),xcart(icart,2)) - sizemax(3) = dmax1(sizemax(3),xcart(icart,3)) - end do - end do - end do - icart = 0 - do itype = 1, ntype - do imol = 1, nmols(itype) - do iatom = 1, natoms(itype) - icart = icart + 1 - sizemin(1) = dmin1(sizemin(1),xcart(icart,1)) - sizemin(2) = dmin1(sizemin(2),xcart(icart,2)) - sizemin(3) = dmin1(sizemin(3),xcart(icart,3)) - sizemax(1) = dmax1(sizemax(1),xcart(icart,1)) - sizemax(2) = dmax1(sizemax(2),xcart(icart,2)) - sizemax(3) = dmax1(sizemax(3),xcart(icart,3)) - end do - end do - end do - write(*,*) ' Mininum and maximum coordinates after constraint fitting: ' - write(*,*) ' x: [ ', sizemin(1),', ', sizemax(1), ' ] ' - write(*,*) ' y: [ ', sizemin(2),', ', sizemax(2), ' ] ' - write(*,*) ' z: [ ', sizemin(3),', ', sizemax(3), ' ] ' - sizemin = sizemin - 1.1d0 * radmax - sizemax = sizemax + 1.1d0 * radmax - ! When *not* using PBC, actually PBCs are used, but particles - ! wont interact across the PBC because there is an empty layer - ! of cells, given by the shift in the sizemin and sizemax dimensions, - ! above defined. - if (.not. using_pbc) then - pbc_min = sizemin - pbc_max = sizemax - pbc_length = pbc_max - pbc_min - end if - - ! Computing the size of the patches - write(*,*) ' Computing size of patches... ' - cell_side = discale * (1.01d0 * radmax) - do i = 1, 3 - ncells(i) = max(1,floor(pbc_length(i)/cell_side)) - cell_length(i) = pbc_length(i) / ncells(i) - end do - write(*,*) ' Number of cells in each direction and cell sides: ' - write(*,*) ' x: ', ncells(1), ' cells of size ', cell_length(1) - write(*,*) ' y: ', ncells(2), ' cells of size ', cell_length(2) - write(*,*) ' z: ', ncells(3), ' cells of size ', cell_length(3) - write(*,"(a, 3(f10.5))") ' Cell-system length: ', pbc_length(1:3) - - ! Allocate arrays depending on the number of cells - allocate(latomfirst(ncells(1),ncells(2),ncells(3))) - allocate(latomfix(ncells(1),ncells(2),ncells(3))) - allocate(lcellnext(ncells(1)*ncells(2)*ncells(3))) - allocate(empty_cell(ncells(1),ncells(2),ncells(3))) - - ! Reseting linked lists arrays - latomfix(:,:,:) = 0 - latomfirst(:,:,:) = 0 - latomnext(:) = 0 - empty_cell(:,:,:) = .true. - - ! If there are fixed molecules, add them permanently to the latomfix array - if(fix) then - write(*,*) ' Add fixed molecules to permanent arrays... ' - icart = ntotat - nfixedat - do iftype = ntype + 1, ntype_with_fixed - do ifatom = 1, natoms(iftype) - icart = icart + 1 - call setcell(xcart(icart,:),cell) - latomnext(icart) = latomfix(cell(1),cell(2),cell(3)) - latomfix(cell(1),cell(2),cell(3)) = icart - latomfirst(cell(1),cell(2),cell(3)) = icart - ibtype(icart) = iftype - ibmol(icart) = 1 - end do - end do - end if - - ! Reseting mass centers to be within the regions - - write(*,*) ' Reseting center of mass... ' - cm_min = 1.d20 - cm_max = -1.d20 - icart = 0 - do itype = 1, ntype - do imol = 1, nmols(itype) - xcm(:) = 0.d0 - do iatom = 1, natoms(itype) - icart = icart + 1 - xcm(:) = xcm(:) + xcart(icart,1:3) - end do - xcm(:) = xcm(:) / natoms(itype) - cm_min(itype,:) = dmin1(cm_min(itype,:),xcm(:)) - cm_max(itype,:) = dmax1(cm_max(itype,:),xcm(:)) - end do - end do - - ! If there is a restart file for all system, read it - - if ( restart_from(0) /= 'none' ) then - record = restart_from(0) - write(*,*) ' Restarting all system from file: ', trim(adjustl(record)) - open(10,file=restart_from(0),status='old',action='read',iostat=ioerr) - ilubar = 0 - ilugan = ntotmol*3 - do i = 1, ntotmol - read(10,*,iostat=ioerr) x(ilubar+1), x(ilubar+2), x(ilubar+3), & - x(ilugan+1), x(ilugan+2), x(ilugan+3) - if ( ioerr /= 0 ) then - write(*,*) ' ERROR: Could not read restart file: ', trim(adjustl(record)) - stop exit_code_open_file - end if - ilubar = ilubar + 3 - ilugan = ilugan + 3 - end do - close(10) - return - end if - - ! Building random initial point - - write(*,dash3_line) - write(*,*) ' Setting initial trial coordinates ... ' - write(*,dash2_line) - - if ( chkgrad ) then - write(*,*) ' For checking gradient, will set avoidoverlap to false. ' - avoidoverlap = .false. - max_guess_try = 1 - scale_rnd_guess = 1.5d0 ! move guess outside the constraints/box - else - max_guess_try = 20 - scale_rnd_guess = 1.d0 - end if - - ! Setting random center of mass coordinates, within size limits - - ilubar = 0 - do itype = 1, ntype - if ( restart_from(itype) /= 'none' ) then - ilubar = ilubar + nmols(itype)*3 - cycle - end if - do imol = 1, nmols(itype) - if ( .not. avoidoverlap ) then - fx = 1.d0 - ntry = 0 - do while((fx > precision) .and. (ntry < max_guess_try)) - ntry = ntry + 1 - call random_number(xrnd) - x(ilubar+1:ilubar+3) = & - cm_min(itype,:) + xrnd(:)*(cm_max(itype,:)-cm_min(itype,:)) * scale_rnd_guess - call restmol(itype,ilubar,n,x,fx,.false.) - end do - else - fx = 1.d0 - ntry = 0 - overlap = .false. - do while((overlap .or. fx > precision) .and. (ntry < max_guess_try)) - overlap = .false. - ntry = ntry + 1 - call random_number(xrnd) - x(ilubar+1:ilubar+3) = & - cm_min(itype,:) + xrnd(:)*(cm_max(itype,:)-cm_min(itype,:)) * scale_rnd_guess - if(fix) then - call setcell(x(ilubar+1:ilubar+3),cell) - icell: do ic = -1, 1 - do jc = -1, 1 - do kc = -1, 1 - if(latomfix(cell_ind(cell(1)+ic,ncells(1)),& - cell_ind(cell(2)+jc,ncells(2)),& - cell_ind(cell(3)+kc,ncells(3)) & - ) /= 0) then - overlap = .true. - exit icell - end if - end do - end do - end do icell - end if - if(.not.overlap) call restmol(itype,ilubar,n,x,fx,.false.) - end do - end if - ilubar = ilubar + 3 - end do - end do - - ! Setting random angles, except if the rotations were constrained - - ilugan = ntotmol*3 - do itype = 1, ntype - if ( restart_from(itype) /= 'none' ) then - ilugan = ilugan + nmols(itype)*3 - cycle - end if - do imol = 1, nmols(itype) - if ( constrain_rot(itype,1) ) then - x(ilugan+1) = ( rot_bound(itype,1,1) - dabs(rot_bound(itype,1,2)) ) + & - 2.d0*rnd()*dabs(rot_bound(itype,1,2)) - else - x(ilugan+1) = twopi*rnd() - end if - if ( constrain_rot(itype,2) ) then - x(ilugan+2) = ( rot_bound(itype,2,1) - dabs(rot_bound(itype,2,2)) ) + & - 2.d0*rnd()*dabs(rot_bound(itype,2,2)) - else - x(ilugan+2) = twopi*rnd() - end if - if ( constrain_rot(itype,3) ) then - x(ilugan+3) = ( rot_bound(itype,3,1) - dabs(rot_bound(itype,3,2)) ) + & - 2.d0*rnd()*dabs(rot_bound(itype,3,2)) - else - x(ilugan+3) = twopi*rnd() - end if - ilugan = ilugan + 3 - end do - end do - - ! Compare analytical and finite-difference gradients - - if(chkgrad) then - cell_side = discale * (1.01d0 * radmax) - do i = 1, 3 - ncells(i) = max(1,floor(pbc_length(i)/cell_side)) - cell_length(i) = pbc_length(i) / ncells(i) - end do - call comparegrad(n,x) - stop - end if - - ! - ! Reading restart files of specific molecule types, if available - ! - - ilubar = 0 - ilugan = ntotmol*3 - do itype = 1, ntype - if ( restart_from(itype) /= 'none' ) then - record = restart_from(itype) - write(*,dash3_line) - write(*,*) ' Molecules of type: ', input_itype(itype) - write(*,*) ' Will restart coordinates from: ', trim(adjustl(record)) - open(10,file=record,status='old',action='read',iostat=ioerr) - if ( ioerr /= 0 ) then - write(*,*) ' ERROR: Could not open restart file: ', trim(adjustl(record)) - stop exit_code_open_file - end if - do i = 1, nmols(itype) - read(10,*,iostat=ioerr) x(ilubar+1), x(ilubar+2), x(ilubar+3), & - x(ilugan+1), x(ilugan+2), x(ilugan+3) - if ( ioerr /= 0 ) then - write(*,*) ' ERROR: Could not read restart file: ', trim(adjustl(record)) - stop exit_code_open_file - end if - ilubar = ilubar + 3 - ilugan = ilugan + 3 - end do - close(10) - call swaptype(n,x,itype,0) ! Initialize swap arrays - call swaptype(n,x,itype,1) ! Set arrays for this type - call computef(n,x,fx) - write(*,*) ' Maximum violation of the restraints: ', frest - write(*,*) ' Maximum violation of minimum atom distances: ', fdist - call swaptype(n,x,itype,3) ! Restore all-molecule arrays - else - ilubar = ilubar + nmols(itype)*3 - ilugan = ilugan + nmols(itype)*3 - end if - end do - - ! Return with current random point (not default) - - if(randini) return - - ! Adjusting current point to fit the constraints - - init1 = .true. - call swaptype(n,x,itype,0) ! Initialize swap arrays - itype = 0 - do while( itype <= ntype ) - itype = itype + 1 - if ( itype == ntype + 1 ) then - call swaptype(n,x,itype,3) ! Restore arrays for all molecules - exit - end if - if ( restart_from(itype) /= 'none' ) cycle - call swaptype(n,x,itype,1) ! Set arrays for this type - write(*,dash3_line) - write(*,*) ' Molecules of type: ', input_itype(itype) - write(*,*) ' Adjusting random positions to fit the constraints. ' - i = 0 - call computef(n,x,fx) - hasbad = .true. - do while( frest > precision .and. i <= nloop0_type(itype)-1 .and. hasbad) - i = i + 1 - write(*,prog1_line) - call pgencan(n,x,fx) - call computef(n,x,fx) - if(frest > precision) then - write(*,"( a,i6,a,i6 )")' Fixing bad orientations ... ', i,' of ', nloop0_type(itype) - movebadprint = .true. - call movebad(n,x,fx,movebadprint) - end if - end do - write(*,*) ' Restraint-only function value: ', fx - write(*,*) ' Maximum violation of the restraints: ', frest - call swaptype(n,x,itype,2) ! Save results for this type - end do - init1 = .false. - write(*,hash3_line) - - ! Deallocate allocated arrays - deallocate(cm_min, cm_max) - - return -end subroutine initial - +! +! Written by Leandro Martínez, 2009-2011. +! Copyright (c) 2009-2018, Leandro Martínez, Jose Mario Martinez, +! Ernesto G. Birgin. +! +! Subroutine initial: Subroutine that reset parameters and +! builds the initial point +! + +subroutine initial(n,x) + + use exit_codes + use sizes + use cell_indexing, only: setcell + use compute_data, only : ntotmol, ntype, nfixedat, ntotat, ncells, nmols, natoms, idfirst, & + ibmol, ibtype, scale, scale2, fdist, frest, sizemin, sizemax, cell_length, radmax, xcart, coor, & + rot_bound, dmax, constrain_rot, comptype, fixedatom, init1, move, latomnext, latomfirst, latomfix, & + lcellfirst, lcellnext, empty_cell + use input, only : randini, ntype_with_fixed, fix, chkgrad, avoidoverlap,& + discale, precision, sidemax, restart_from, input_itype,& + nloop0_type + use usegencan + use ahestetic + use pbc + implicit none + integer :: n, i, j, idatom, iatom, ilubar, ilugan, icart, itype, & + imol, ntry, cell(3), ic, jc, kc, ifatom, & + idfatom, iftype, jatom, ioerr, max_guess_try + + double precision :: x(n), beta, gamma, theta, & + fx, cell_side, rnd, xrnd(3), v1(3), v2(3), v3(3), xcm(3) + double precision :: scale_rnd_guess + double precision, parameter :: twopi = 8.d0*datan(1.d0) + double precision, allocatable :: cm_min(:,:), cm_max(:,:) + + logical :: overlap, movebadprint, hasbad + + character(len=strl) :: record + + ! We need to initialize the move logical variable + move = .false. + + ! Default status of the function evaluation + init1 = .false. + lcellfirst = 0 + + ! Initialize the comptype logical array + comptype(1:ntype_with_fixed) = .true. + + ! Penalty factors for the objective function relative to restrictions + ! Default values: scale = 1.d2, scale2 = 1.d1 + + scale = 1.d0 + scale2 = 1.d-2 + call tobar() + + ! Compute maximum internal distance within each type of molecule + + do itype = 1, ntype + dmax(itype) = 0.d0 + idatom = idfirst(itype) - 1 + do iatom = 1, natoms(itype) - 1 + do jatom = iatom + 1, natoms(itype) + dmax(itype) = dmax1 ( dmax(itype),& + (coor(idatom+iatom,1)-coor(idatom+jatom,1))**2+& + (coor(idatom+iatom,2)-coor(idatom+jatom,2))**2+& + (coor(idatom+iatom,3)-coor(idatom+jatom,3))**2 ) + end do + end do + dmax(itype) = dsqrt(dmax(itype)) + write(*,*) ' Maximum internal distance of type ',itype,': ',& + dmax(itype) + if(dmax(itype).eq.0.) dmax(itype) = 1.d0 + end do + + ! Maximum size of the system: if you system is very large (about + ! 80 nm wide), increase the sidemax parameter. + ! Otherwise, the packing can be slow and unsucesful + allocate(cm_min(ntype,3), cm_max(ntype,3)) + x(:) = 0.d0 + call restmol(1,0,n,x,fx,.true.) + sizemin(1:3) = x(1:3) - sidemax + sizemax(1:3) = x(1:3) + sidemax + write(*,*) ' All atoms must be within these coordinates: ' + write(*,*) ' x: [ ', sizemin(1),', ', sizemax(1), ' ] ' + write(*,*) ' y: [ ', sizemin(2),', ', sizemax(2), ' ] ' + write(*,*) ' z: [ ', sizemin(3),', ', sizemax(3), ' ] ' + write(*,*) ' If the system is larger than this, increase the sidemax parameter. ' + + ! Create first aleatory guess + + i = 0 + j = ntotmol*3 + do itype = 1, ntype + do imol = 1, nmols(itype) + x(i+1) = sizemin(1) + rnd()*(sizemax(1)-sizemin(1)) + x(i+2) = sizemin(2) + rnd()*(sizemax(2)-sizemin(2)) + x(i+3) = sizemin(3) + rnd()*(sizemax(3)-sizemin(3)) + if ( constrain_rot(itype,1) ) then + x(j+1) = ( rot_bound(itype,1,1) - dabs(rot_bound(itype,1,2)) ) + & + 2.d0*rnd()*dabs(rot_bound(itype,1,2)) + else + x(j+1) = twopi*rnd() + end if + if ( constrain_rot(itype,2) ) then + x(j+2) = ( rot_bound(itype,2,1) - dabs(rot_bound(itype,2,2)) ) + & + 2.d0*rnd()*dabs(rot_bound(itype,2,2)) + else + x(j+2) = twopi*rnd() + end if + if ( constrain_rot(itype,3) ) then + x(j+3) = ( rot_bound(itype,3,1) - dabs(rot_bound(itype,3,2)) ) + & + 2.d0*rnd()*dabs(rot_bound(itype,3,2)) + else + x(j+3) = twopi*rnd() + end if + i = i + 3 + j = j + 3 + end do + end do + + ! Initialize cartesian coordinate array for the first time + + ilubar = 0 + ilugan = ntotmol*3 + icart = 0 + do itype = 1, ntype + do imol = 1, nmols(itype) + xcm = x(ilubar+1:ilubar+3) + beta = x(ilugan+1) + gamma = x(ilugan+2) + theta = x(ilugan+3) + call eulerrmat(beta,gamma,theta,v1,v2,v3) + idatom = idfirst(itype) - 1 + do iatom = 1, natoms(itype) + icart = icart + 1 + idatom = idatom + 1 + call compcart(xcart(icart,1:3),xcm,coor(idatom,1:3),v1,v2,v3) + fixedatom(icart) = .false. + + end do + end do + end do + if(fix) then + icart = ntotat - nfixedat + do iftype = ntype + 1, ntype_with_fixed + idfatom = idfirst(iftype) - 1 + do ifatom = 1, natoms(iftype) + idfatom = idfatom + 1 + icart = icart + 1 + xcart(icart,1) = coor(idfatom,1) + xcart(icart,2) = coor(idfatom,2) + xcart(icart,3) = coor(idfatom,3) + fixedatom(icart) = .true. + + ! Check if fixed molecules are compatible with PBC given + if (using_pbc) then + do i = 1, 3 + if (xcart(icart, i) < pbc_min(i) .or. xcart(icart, i) > pbc_max(i)) then + write(*,*) "ERROR: Fixed molecule are outside the PBC box:" + write(*,*) " Atom: ", ifatom, " of molecule: ", input_itype(iftype), " - coordinate: ", i + write(*,*) " ", xcart(icart, i), " not in [", pbc_min(i), ", ", pbc_max(i), "]" + write(*,*) "(after translating/rotation the fixed molecule with the given orientation)" + stop exit_code_input_error + end if + end do + end if + end do + end do + end if + + ! Performing some steps of optimization for the restrictions only + + write(*,hash3_line) + write(*,"(' Building initial approximation ... ' )") + write(*,hash3_line) + write(*,"(' Adjusting initial point to fit the constraints ')") + write(*,dash2_line) + init1 = .true. + call swaptype(n,x,itype,0) ! Initialize swap arrays + itype = 0 + do while( itype <= ntype ) + itype = itype + 1 + if ( itype <= ntype ) then + call swaptype(n,x,itype,1) ! Set arrays for this type + else + call swaptype(n,x,itype,3) ! Restore arrays if itype = ntype + 1 + exit + end if + write(*,dash3_line) + write(*,*) ' Molecules of type: ', input_itype(itype) + write(*,*) + i = 0 + hasbad = .true. + call computef(n,x,fx) + do while( frest > precision .and. i.le. nloop0_type(itype)-1 .and. hasbad) + i = i + 1 + write(*,prog1_line) + call pgencan(n,x,fx) + call computef(n,x,fx) + if(frest > precision) then + write(*,"( a,i6,a,i6 )")' Fixing bad orientations ... ', i,' of ', nloop0_type(itype) + movebadprint = .true. + call movebad(n,x,fx,movebadprint) + end if + end do + write(*,*) + write(*,*) ' Restraint-only function value: ', fx + write(*,*) ' Maximum violation of the restraints: ', frest + + call swaptype(n,x,itype,2) ! Save current type results + + if( hasbad .and. frest > precision ) then + write(*,*) ' ERROR: Packmol was unable to put the molecules' + write(*,*) ' in the desired regions even without' + write(*,*) ' considering distance tolerances. ' + write(*,*) ' Probably there is something wrong with' + write(*,*) ' the constraints, since it seems that' + write(*,*) ' the molecules cannot satisfy them at' + write(*,*) ' at all. ' + write(*,*) ' Please check the spatial constraints and' + write(*,*) ' try again.' + if ( i .ge. nloop0_type(itype)-1 ) then + end if + write(*,*) ' >The maximum number of cycles (',nloop0_type(itype),') was achieved.' + write(*,*) ' You may try increasing it with the',' nloop0 keyword, as in: nloop0 1000 ' + stop exit_code_failed_to_converge + end if + end do + init1 = .false. + + ! Rescaling sizemin and sizemax in order to build the patch of cells + + write(*,dash3_line) + write(*,*) ' Rescaling maximum and minimum coordinates... ' + sizemin(1:3) = 1.d20 + sizemax(1:3) = -1.d20 + ! Maximum and minimum coordinates of fixed molecules + icart = ntotat - nfixedat + do itype = ntype + 1, ntype_with_fixed + do imol = 1, nmols(itype) + do iatom = 1, natoms(itype) + icart = icart + 1 + sizemin(1) = dmin1(sizemin(1),xcart(icart,1)) + sizemin(2) = dmin1(sizemin(2),xcart(icart,2)) + sizemin(3) = dmin1(sizemin(3),xcart(icart,3)) + sizemax(1) = dmax1(sizemax(1),xcart(icart,1)) + sizemax(2) = dmax1(sizemax(2),xcart(icart,2)) + sizemax(3) = dmax1(sizemax(3),xcart(icart,3)) + end do + end do + end do + icart = 0 + do itype = 1, ntype + do imol = 1, nmols(itype) + do iatom = 1, natoms(itype) + icart = icart + 1 + sizemin(1) = dmin1(sizemin(1),xcart(icart,1)) + sizemin(2) = dmin1(sizemin(2),xcart(icart,2)) + sizemin(3) = dmin1(sizemin(3),xcart(icart,3)) + sizemax(1) = dmax1(sizemax(1),xcart(icart,1)) + sizemax(2) = dmax1(sizemax(2),xcart(icart,2)) + sizemax(3) = dmax1(sizemax(3),xcart(icart,3)) + end do + end do + end do + write(*,*) ' Mininum and maximum coordinates after constraint fitting: ' + write(*,*) ' x: [ ', sizemin(1),', ', sizemax(1), ' ] ' + write(*,*) ' y: [ ', sizemin(2),', ', sizemax(2), ' ] ' + write(*,*) ' z: [ ', sizemin(3),', ', sizemax(3), ' ] ' + sizemin = sizemin - 1.1d0 * radmax + sizemax = sizemax + 1.1d0 * radmax + ! When *not* using PBC, actually PBCs are used, but particles + ! wont interact across the PBC because there is an empty layer + ! of cells, given by the shift in the sizemin and sizemax dimensions, + ! above defined. + if (.not. using_pbc) then + pbc_min = sizemin + pbc_max = sizemax + pbc_length = pbc_max - pbc_min + end if + + ! Computing the size of the patches + write(*,*) ' Computing size of patches... ' + cell_side = discale * (1.01d0 * radmax) + do i = 1, 3 + ncells(i) = max(1,floor(pbc_length(i)/cell_side)) + cell_length(i) = pbc_length(i) / ncells(i) + end do + write(*,*) ' Number of cells in each direction and cell sides: ' + write(*,*) ' x: ', ncells(1), ' cells of size ', cell_length(1) + write(*,*) ' y: ', ncells(2), ' cells of size ', cell_length(2) + write(*,*) ' z: ', ncells(3), ' cells of size ', cell_length(3) + write(*,"(a, 3(f10.5))") ' Cell-system length: ', pbc_length(1:3) + + ! Allocate arrays depending on the number of cells + allocate(latomfirst(ncells(1),ncells(2),ncells(3))) + allocate(latomfix(ncells(1),ncells(2),ncells(3))) + allocate(lcellnext(ncells(1)*ncells(2)*ncells(3))) + allocate(empty_cell(ncells(1),ncells(2),ncells(3))) + + ! Reseting linked lists arrays + latomfix(:,:,:) = 0 + latomfirst(:,:,:) = 0 + latomnext(:) = 0 + empty_cell(:,:,:) = .true. + + ! If there are fixed molecules, add them permanently to the latomfix array + if(fix) then + write(*,*) ' Add fixed molecules to permanent arrays... ' + icart = ntotat - nfixedat + do iftype = ntype + 1, ntype_with_fixed + do ifatom = 1, natoms(iftype) + icart = icart + 1 + call setcell(xcart(icart,:),cell) + latomnext(icart) = latomfix(cell(1),cell(2),cell(3)) + latomfix(cell(1),cell(2),cell(3)) = icart + latomfirst(cell(1),cell(2),cell(3)) = icart + ibtype(icart) = iftype + ibmol(icart) = 1 + end do + end do + end if + + ! Reseting mass centers to be within the regions + + write(*,*) ' Reseting center of mass... ' + cm_min = 1.d20 + cm_max = -1.d20 + icart = 0 + do itype = 1, ntype + do imol = 1, nmols(itype) + xcm(:) = 0.d0 + do iatom = 1, natoms(itype) + icart = icart + 1 + xcm(:) = xcm(:) + xcart(icart,1:3) + end do + xcm(:) = xcm(:) / natoms(itype) + cm_min(itype,:) = dmin1(cm_min(itype,:),xcm(:)) + cm_max(itype,:) = dmax1(cm_max(itype,:),xcm(:)) + end do + end do + + ! If there is a restart file for all system, read it + + if ( restart_from(0) /= 'none' ) then + record = restart_from(0) + write(*,*) ' Restarting all system from file: ', trim(adjustl(record)) + open(10,file=restart_from(0),status='old',action='read',iostat=ioerr) + ilubar = 0 + ilugan = ntotmol*3 + do i = 1, ntotmol + read(10,*,iostat=ioerr) x(ilubar+1), x(ilubar+2), x(ilubar+3), & + x(ilugan+1), x(ilugan+2), x(ilugan+3) + if ( ioerr /= 0 ) then + write(*,*) ' ERROR: Could not read restart file: ', trim(adjustl(record)) + stop exit_code_open_file + end if + ilubar = ilubar + 3 + ilugan = ilugan + 3 + end do + close(10) + return + end if + + ! Building random initial point + + write(*,dash3_line) + write(*,*) ' Setting initial trial coordinates ... ' + write(*,dash2_line) + + if ( chkgrad ) then + write(*,*) ' For checking gradient, will set avoidoverlap to false. ' + avoidoverlap = .false. + max_guess_try = 1 + scale_rnd_guess = 1.5d0 ! move guess outside the constraints/box + else + max_guess_try = 20 + scale_rnd_guess = 1.d0 + end if + + ! Setting random center of mass coordinates, within size limits + + ilubar = 0 + do itype = 1, ntype + if ( restart_from(itype) /= 'none' ) then + ilubar = ilubar + nmols(itype)*3 + cycle + end if + do imol = 1, nmols(itype) + if ( .not. avoidoverlap ) then + fx = 1.d0 + ntry = 0 + do while((fx > precision) .and. (ntry < max_guess_try)) + ntry = ntry + 1 + call random_number(xrnd) + x(ilubar+1:ilubar+3) = & + cm_min(itype,:) + xrnd(:)*(cm_max(itype,:)-cm_min(itype,:)) * scale_rnd_guess + call restmol(itype,ilubar,n,x,fx,.false.) + end do + else + fx = 1.d0 + ntry = 0 + overlap = .false. + do while((overlap .or. fx > precision) .and. (ntry < max_guess_try)) + overlap = .false. + ntry = ntry + 1 + call random_number(xrnd) + x(ilubar+1:ilubar+3) = & + cm_min(itype,:) + xrnd(:)*(cm_max(itype,:)-cm_min(itype,:)) * scale_rnd_guess + if(fix) then + call setcell(x(ilubar+1:ilubar+3),cell) + icell: do ic = -1, 1 + do jc = -1, 1 + do kc = -1, 1 + if(latomfix(cell_ind(cell(1)+ic,ncells(1)),& + cell_ind(cell(2)+jc,ncells(2)),& + cell_ind(cell(3)+kc,ncells(3)) & + ) /= 0) then + overlap = .true. + exit icell + end if + end do + end do + end do icell + end if + if(.not.overlap) call restmol(itype,ilubar,n,x,fx,.false.) + end do + end if + ilubar = ilubar + 3 + end do + end do + + ! Setting random angles, except if the rotations were constrained + + ilugan = ntotmol*3 + do itype = 1, ntype + if ( restart_from(itype) /= 'none' ) then + ilugan = ilugan + nmols(itype)*3 + cycle + end if + do imol = 1, nmols(itype) + if ( constrain_rot(itype,1) ) then + x(ilugan+1) = ( rot_bound(itype,1,1) - dabs(rot_bound(itype,1,2)) ) + & + 2.d0*rnd()*dabs(rot_bound(itype,1,2)) + else + x(ilugan+1) = twopi*rnd() + end if + if ( constrain_rot(itype,2) ) then + x(ilugan+2) = ( rot_bound(itype,2,1) - dabs(rot_bound(itype,2,2)) ) + & + 2.d0*rnd()*dabs(rot_bound(itype,2,2)) + else + x(ilugan+2) = twopi*rnd() + end if + if ( constrain_rot(itype,3) ) then + x(ilugan+3) = ( rot_bound(itype,3,1) - dabs(rot_bound(itype,3,2)) ) + & + 2.d0*rnd()*dabs(rot_bound(itype,3,2)) + else + x(ilugan+3) = twopi*rnd() + end if + ilugan = ilugan + 3 + end do + end do + + ! Compare analytical and finite-difference gradients + + if(chkgrad) then + cell_side = discale * (1.01d0 * radmax) + do i = 1, 3 + ncells(i) = max(1,floor(pbc_length(i)/cell_side)) + cell_length(i) = pbc_length(i) / ncells(i) + end do + call comparegrad(n,x) + stop + end if + + ! + ! Reading restart files of specific molecule types, if available + ! + + ilubar = 0 + ilugan = ntotmol*3 + do itype = 1, ntype + if ( restart_from(itype) /= 'none' ) then + record = restart_from(itype) + write(*,dash3_line) + write(*,*) ' Molecules of type: ', input_itype(itype) + write(*,*) ' Will restart coordinates from: ', trim(adjustl(record)) + open(10,file=record,status='old',action='read',iostat=ioerr) + if ( ioerr /= 0 ) then + write(*,*) ' ERROR: Could not open restart file: ', trim(adjustl(record)) + stop exit_code_open_file + end if + do i = 1, nmols(itype) + read(10,*,iostat=ioerr) x(ilubar+1), x(ilubar+2), x(ilubar+3), & + x(ilugan+1), x(ilugan+2), x(ilugan+3) + if ( ioerr /= 0 ) then + write(*,*) ' ERROR: Could not read restart file: ', trim(adjustl(record)) + stop exit_code_open_file + end if + ilubar = ilubar + 3 + ilugan = ilugan + 3 + end do + close(10) + call swaptype(n,x,itype,0) ! Initialize swap arrays + call swaptype(n,x,itype,1) ! Set arrays for this type + call computef(n,x,fx) + write(*,*) ' Maximum violation of the restraints: ', frest + write(*,*) ' Maximum violation of minimum atom distances: ', fdist + call swaptype(n,x,itype,3) ! Restore all-molecule arrays + else + ilubar = ilubar + nmols(itype)*3 + ilugan = ilugan + nmols(itype)*3 + end if + end do + + ! Return with current random point (not default) + + if(randini) return + + ! Adjusting current point to fit the constraints + + init1 = .true. + call swaptype(n,x,itype,0) ! Initialize swap arrays + itype = 0 + do while( itype <= ntype ) + itype = itype + 1 + if ( itype == ntype + 1 ) then + call swaptype(n,x,itype,3) ! Restore arrays for all molecules + exit + end if + if ( restart_from(itype) /= 'none' ) cycle + call swaptype(n,x,itype,1) ! Set arrays for this type + write(*,dash3_line) + write(*,*) ' Molecules of type: ', input_itype(itype) + write(*,*) ' Adjusting random positions to fit the constraints. ' + i = 0 + call computef(n,x,fx) + hasbad = .true. + do while( frest > precision .and. i <= nloop0_type(itype)-1 .and. hasbad) + i = i + 1 + write(*,prog1_line) + call pgencan(n,x,fx) + call computef(n,x,fx) + if(frest > precision) then + write(*,"( a,i6,a,i6 )")' Fixing bad orientations ... ', i,' of ', nloop0_type(itype) + movebadprint = .true. + call movebad(n,x,fx,movebadprint) + end if + end do + write(*,*) ' Restraint-only function value: ', fx + write(*,*) ' Maximum violation of the restraints: ', frest + call swaptype(n,x,itype,2) ! Save results for this type + end do + init1 = .false. + write(*,hash3_line) + + ! Deallocate allocated arrays + deallocate(cm_min, cm_max) + + return +end subroutine initial + diff --git a/src/output.f90 b/src/output.f90 index 3e530fb..c061a1e 100644 --- a/src/output.f90 +++ b/src/output.f90 @@ -10,7 +10,7 @@ subroutine output(n, x, output_file_name) use exit_codes use sizes - use compute_data + use compute_data, only : ntotmol, ntype, ntotat, nmols, natoms, idfirst, xcart, coor use input use pbc @@ -21,6 +21,7 @@ subroutine output(n, x, output_file_name) icart, i_ref_atom, ioerr, ifirst_mol integer :: nres, imark integer :: i_fixed, i_not_fixed + logical :: has_residue_records double precision :: x(n) double precision :: v1(3), v2(3), v3(3) @@ -328,13 +329,15 @@ subroutine output(n, x, output_file_name) open(15,file=pdbfile(i_fixed),status='old') ifres = 0 + ilres = 0 + has_residue_records = .false. do read(15,str_format,iostat=ioerr) record if ( ioerr /= 0 ) exit if ( record(1:4).eq.'ATOM'.or.record(1:6).eq.'HETATM' ) then read(record(23:26),*,iostat=ioerr) imark if ( ioerr /= 0 ) then - record = pdbfile(i_not_fixed) + record = pdbfile(i_fixed) write(*,*) ' ERROR: Failed reading residue number ',& ' from PDB file: ', trim(adjustl(record)) write(*,*) ' Residue numbers are integers that must',& @@ -348,8 +351,17 @@ subroutine output(n, x, output_file_name) end if if ( ifres .eq. 0 ) ifres = imark ilres = imark + has_residue_records = .true. end if end do + if ( .not. has_residue_records ) then + record = pdbfile(i_fixed) + write(*,*) ' ERROR: Could not determine residue numbering '& + //'for fixed molecule file: ', trim(adjustl(record)) + write(*,*) ' At least one ATOM or HETATM record must be '& + //'present in the input PDB file. ' + stop exit_code_input_error + end if nres = ilres - ifres + 1 iimol = iimol + 1 diff --git a/src/pgencan.f90 b/src/pgencan.f90 index f907298..7e35139 100644 --- a/src/pgencan.f90 +++ b/src/pgencan.f90 @@ -11,8 +11,8 @@ subroutine pgencan(n,x,fx) use sizes - use compute_data - use usegencan + use compute_data, only : ntype, nmols, constrain_rot, rot_bound, init1 + use usegencan, only : maxit, iprint1, iprint2, wi, wd, l, u, g implicit none double precision :: lambda(1), rho(1) diff --git a/src/resetcells.f90 b/src/resetcells.f90 index f3e9ba7..9674be6 100644 --- a/src/resetcells.f90 +++ b/src/resetcells.f90 @@ -1,43 +1,51 @@ -! -! Written by Leandro Martínez, 2009-2011. -! Copyright (c) 2009-2018, Leandro Martínez, Jose Mario Martinez, -! Ernesto G. Birgin. -! -! Subroutine resetcells: Subroutine that resets the occupancy of -! linked cell cells -! - -subroutine resetcells() - use input, only: fix, ntype_with_fixed - use cell_indexing, only: icell_to_cell, index_cell, setcell - use compute_data, only: latomfirst, lcellfirst, lcellnext, empty_cell, ncells, & - ntype, natoms, ntotat, nfixedat, xcart, latomnext - implicit none - integer :: cell(3), icell, iftype, icart, ifatom - - ! Reset data for cells that contain fixed atom - lcellnext(:) = 0 - lcellfirst = 0 - latomfirst(:,:,:) = 0 - latomnext(:) = 0 - empty_cell(:,:,:) = .true. - if(fix) then - icart = ntotat - nfixedat - do iftype = ntype + 1, ntype_with_fixed - do ifatom = 1, natoms(iftype) - icart = icart + 1 - call setcell(xcart(icart,:),cell) - latomnext(icart) = latomfirst(cell(1),cell(2),cell(3)) - latomfirst(cell(1),cell(2),cell(3)) = icart - if ( empty_cell(cell(1),cell(2),cell(3)) ) then - empty_cell(cell(1),cell(2),cell(3)) = .false. - icell = index_cell(cell,ncells) - lcellnext(icell) = lcellfirst - lcellfirst = icell - end if - end do - end do - end if - -end subroutine resetcells - +! +! Written by Leandro Martínez, 2009-2011. +! Copyright (c) 2009-2018, Leandro Martínez, Jose Mario Martinez, +! Ernesto G. Birgin. +! +! Subroutine resetcells: Subroutine that resets the occupancy of +! linked cell cells +! + +subroutine resetcells() + use input, only: fix, ntype_with_fixed + use cell_indexing, only: icell_to_cell, index_cell, setcell + use compute_data, only: latomfirst, lcellfirst, lcellnext, empty_cell, ncells, & + ntype, natoms, ntotat, nfixedat, xcart, latomnext + implicit none + integer :: cell(3), icell, next_icell, iftype, icart, ifatom + + ! Incrementally reset only previously-occupied cells instead of + ! clearing all ncells^3 entries. Walk the old occupied-cell linked + ! list and reset each cell's data, then rebuild the fixed-atom lists. + icell = lcellfirst + do while ( icell > 0 ) + next_icell = lcellnext(icell) + call icell_to_cell(icell, ncells, cell) + latomfirst(cell(1),cell(2),cell(3)) = 0 + empty_cell(cell(1),cell(2),cell(3)) = .true. + lcellnext(icell) = 0 + icell = next_icell + end do + lcellfirst = 0 + + if(fix) then + icart = ntotat - nfixedat + do iftype = ntype + 1, ntype_with_fixed + do ifatom = 1, natoms(iftype) + icart = icart + 1 + call setcell(xcart(icart,:),cell) + latomnext(icart) = latomfirst(cell(1),cell(2),cell(3)) + latomfirst(cell(1),cell(2),cell(3)) = icart + + if ( empty_cell(cell(1),cell(2),cell(3)) ) then + empty_cell(cell(1),cell(2),cell(3)) = .false. + icell = index_cell(cell,ncells) + lcellnext(icell) = lcellfirst + lcellfirst = icell + end if + end do + end do + end if + +end subroutine resetcells diff --git a/src/restmol.f90 b/src/restmol.f90 index 909a41c..6e6d5ee 100644 --- a/src/restmol.f90 +++ b/src/restmol.f90 @@ -12,7 +12,7 @@ subroutine restmol(itype,ilubar,n,x,fx,solve) use sizes - use compute_data + use compute_data, only : ntotmol, ntype, nmols, comptype, init1, xmol, compsafe use usegencan implicit none diff --git a/src/setsizes.f90 b/src/setsizes.f90 index a73337e..ba54b34 100644 --- a/src/setsizes.f90 +++ b/src/setsizes.f90 @@ -325,7 +325,8 @@ subroutine setsizes() ! Allocate arrays depending on ntotat, nn, maxrest, and mrperatom allocate(nratom(ntotat),iratom(ntotat,mrperatom),ibmol(ntotat),& - ibtype(ntotat),xcart(ntotat,3),coor(ntotat,3),& + ibtype(ntotat),fixed_short_marker(ntotat),& + xcart(ntotat,3),coor(ntotat,3),& radius(ntotat),radius_ini(ntotat),fscale(ntotat),& use_short_radius(ntotat), short_radius(ntotat), short_radius_scale(ntotat),& gxcar(ntotat,3),& diff --git a/testing/check_numerics_profiles.sh b/testing/check_numerics_profiles.sh new file mode 100755 index 0000000..09edf70 --- /dev/null +++ b/testing/check_numerics_profiles.sh @@ -0,0 +1,90 @@ +#!/usr/bin/env bash +set -euo pipefail + +ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" +TEST_DIR="${ROOT_DIR}/testing" +TMP_DIR="${TEST_DIR}/tmp_profile_check" + +BASELINE_BIN="${ROOT_DIR}/packmol-baseline" +PERF_BIN="${ROOT_DIR}/packmol-perf-native" +INPUT_FILE="${TEST_DIR}/input_files/benzene2.inp" +INPUT_CHKGRAD="${TMP_DIR}/benzene2_chkgrad.inp" + +OBJ_TOL="${OBJ_TOL:-1e-8}" +GRAD_TOL="${GRAD_TOL:-1e-3}" + +extract_function_value() { + awk '/Function Value =/ {print $4}' "$1" +} + +extract_max_grad_error() { + awk '/Maximum difference =/ {print $7}' "$1" +} + +run_profile() { + local profile="$1" + local binary="$2" + local outdir="${TMP_DIR}/${profile}" + + rm -rf "${outdir}" + mkdir -p "${outdir}" + cp -r "${TEST_DIR}/structure_files" "${outdir}/" + + ( + cd "${ROOT_DIR}" + make clean >/dev/null + make "${profile}" >/dev/null + cp packmol "${binary}" + ) + + ( + cd "${outdir}" + "${binary}" < "${INPUT_CHKGRAD}" > profile.log + ) + + cp "${outdir}/chkgrad.log" "${TMP_DIR}/chkgrad-${profile}.log" +} + +mkdir -p "${TMP_DIR}" +cp "${INPUT_FILE}" "${INPUT_CHKGRAD}" +echo "chkgrad" >> "${INPUT_CHKGRAD}" + +echo "[1/2] Running baseline profile" +run_profile baseline "${BASELINE_BIN}" + +echo "[2/2] Running perf-native profile" +run_profile perf-native "${PERF_BIN}" + +base_obj="$(extract_function_value "${TMP_DIR}/chkgrad-baseline.log")" +perf_obj="$(extract_function_value "${TMP_DIR}/chkgrad-perf-native.log")" +base_grad_err="$(extract_max_grad_error "${TMP_DIR}/chkgrad-baseline.log")" +perf_grad_err="$(extract_max_grad_error "${TMP_DIR}/chkgrad-perf-native.log")" + +python - "$base_obj" "$perf_obj" "$base_grad_err" "$perf_grad_err" "$OBJ_TOL" "$GRAD_TOL" <<'PY' +import math +import sys + +base_obj = float(sys.argv[1]) +perf_obj = float(sys.argv[2]) +base_grad_err = float(sys.argv[3]) +perf_grad_err = float(sys.argv[4]) +obj_tol = float(sys.argv[5]) +grad_tol = float(sys.argv[6]) + +obj_delta = abs(base_obj - perf_obj) +grad_delta = abs(base_grad_err - perf_grad_err) + +print(f"Baseline objective: {base_obj:.16e}") +print(f"Perf-native objective: {perf_obj:.16e}") +print(f"|Δ objective| = {obj_delta:.6e} (tol={obj_tol:.6e})") +print(f"Baseline max grad check error: {base_grad_err:.6e}") +print(f"Perf-native max grad check error: {perf_grad_err:.6e}") +print(f"|Δ grad error| = {grad_delta:.6e} (tol={grad_tol:.6e})") + +if obj_delta > obj_tol: + raise SystemExit("Objective delta exceeds tolerance") +if grad_delta > grad_tol: + raise SystemExit("Gradient delta exceeds tolerance") +PY + +echo "Numerics comparison passed."