1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8 module neighbour 9 use precision, only: dp 10 use sys, only: die 11 use alloc, only: re_alloc, de_alloc 12 implicit none 13 14 private 15 16 public :: mneighb, reset_neighbour_arrays 17 18 character(len=*), parameter :: myName = 'neighbour ' 19 integer, public :: maxnna = 200 20 integer, pointer, public :: jan(:) 21 real(dp), pointer, public :: r2ij(:) 22 real(dp), pointer, public :: xij(:,:) 23 logical :: pointers_allocated = .false. 24 integer :: maxna = -1 25 integer :: maxnm = -1 26 integer :: maxnnm = -1 27 integer :: maxnem = -1 28 integer, parameter :: nx = 3 ! Define space dimension 29 30 integer :: INX(NX), I1NX(NX), I2NX(NX), 31 & J1NX(NX), J2NX(NX), I1EMX(NX), 32 & I2EMX(NX), IMX(NX), I1MX(NX), 33 & I2MX(NX), NEMX(NX), NMX(NX), NNX(NX) 34 real(dp) :: DMX(NX), DX(NX), DX0M(NX), 35 & CELMSH(NX*NX), RCELL(NX*NX), 36 & RMCELL(NX*NX) 37 integer, pointer :: IANEXT(:), IAPREV(:), IEMA(:), 38 & IA1M(:), IDNM(:), IMESH(:) 39 real(dp), pointer :: DXAM(:,:), DXNM(:,:) 40 real(dp) :: celast(nx,nx) = 0.0_dp, 41 & rglast = 0.0_dp 42 real(dp), public :: x0(nx) 43 44 contains 45 46 subroutine mneighb( cell, range, na, xa, ia, isc, 47 & nna ) 48 49C ******************************************************************** 50C Finds the neighbours of an atom in a cell with periodic boundary 51C conditions. This is an interface to routine ranger, which has 52C extended functionalities. 53C Written by J.M.Soler. March 1997. 54C *********** INPUT ************************************************** 55C REAL*8 CELL(3,3) : Unit cell vectors CELL(IXYZ,IVECT) 56C REAL*8 range : Maximum distance of neighbours required 57C INTEGER NA : Number of atoms 58C REAL*8 XA(3,NA) : Atomic positions in cartesian coordinates 59C INTEGER IA : Atom whose neighbours are needed. 60C A routine initialization must be done by 61C a first call with IA = 0 62C INTEGER ISC : Single-counting switch (0=No, 1=Yes). If ISC=1, 63C only neighbours with JA.LE.IA are included in JAN 64!!!!! C INTEGER MAXNNA : Size of arrays JAN, XIJ and R2IJ 65C *********** OUTPUT ************************************************* 66C INTEGER NNA : Number of neighbour atoms within range of IA 67C INTEGER JAN(NNA) : Atom-index of neighbours 68C REAL*8 XIJ(3,NNA): Vectors from atom IA to neighbours 69C REAL*8 R2IJ(NNA) : Squared distances to neighbours 70C *********** UNITS ************************************************** 71C Units of CELL, range and XA are arbitrary but must be equal 72C *********** SUBROUTINES USED *************************************** 73C CHKDIM, DISMIN, RECLAT 74C *********** BEHAVIOUR ********************************************** 75C CPU time and memory scale linearly with the number of atoms, for 76C sufficiently large numbers. 77C If internal dimension variables are too small, an include file named 78C neighb.h is printed with the required dimensions and the routine 79C stops, asking to be recompiled. Then, neighb.h is automatically 80C included in the new compilation, if the source file is in the same 81C directory where the program has run. Initially, you can make all 82C the parameters in neighb.h equal to 1. 83C Different ranges can be used for different atoms, but for good 84C performance, the largest range should be used in the initial 85C call (with IA=0). 86C There are no limitations regarding cell shape or size. The range may 87C be larger than the cell size, in which case many 'images' of the 88C same atom will be included in the neighbour list, with different 89C interatomic vectors XIJ and distances R2IJ. 90C The atom IA itself is included in the neighbour list, with zero 91C distance. You have to discard it if you want so. 92C If the number of neighbour atoms found is larger than the size of 93C the arrays JAN, XIJ and R2IJ, i.e. if NNAout > NNAin, these arrays 94C are filled only up to their size NNAin. With dynamic memory 95C allocation, this allows to find first the required array sizes 96C and then find the neighbours. Notice however that no warning is 97C given, so that you should always check that NNAout.LE.NNAin. 98C *********** USAGE ************************************************** 99C Sample usage for a molecular dynamics simulation: 100C DIMENSION JAN(MAXNNA), XIJ(3,MAXNNA) 101C Define CELL and initial positions XA 102C DO ITER = 1,NITER (Molecular dynamics iteration) 103C NNA = MAXNNA 104C CALL NEIGHB( CELL, range, NA, XA, 0, 1, NNA, JAN, XIJ, R2IJ ) 105C IF (NNA .GT. MAXNNA) STOP 'MAXNNA too small' 106C Initialize to zero all atomic forces FA(IX,IA) 107C DO IA = 1,NA (Loop on atoms) 108C NNA = MAXNNA 109C CALL NEIGHB( CELL, range, NA, XA, IA, 1, NNA, JAN, XIJ, R2IJ ) 110C IF (NNA .GT. MAXNNA) STOP 'MAXNNA too small' 111C DO IN = 1,NNA (Loop on neighbours of IA) 112C JA = JAN(IN) (Atomic index of neighbour) 113C RIJ = SQRT(R2IJ(IN)) (Interatomic distance) 114C IF (RIJ .GT. 1.D-12) THEN (Discard atom itself) 115C Find interatomic force FIJ( RIJ ) 116C DO IX = 1,3 117C FA(IX,IA) = FA(IX,IA) - FIJ * XIJ(IX,IN) / RIJ 118C FA(IX,JA) = FA(IX,JA) + FIJ * XIJ(IX,IN) / RIJ 119C ENDDO 120C ENDIF 121C ENDDO 122C ENDDO 123C Move atomic positions XA (Molecular dynamics step) 124C ENDDO 125C ******************************************************************** 126 implicit none 127C Argument types and dimensions 128 integer, intent(in) :: ia, isc, na 129 integer, intent(out) :: nna 130 real(dp) :: cell(nx,nx), range, xa(nx,na) 131 132C Internal variables 133 integer, save :: iamove(1) = 0 134 logical, save :: first_time = .true. 135 logical :: samcel 136 integer :: IX, JX 137 138 call sizeup_neighbour_arrays( maxnna ) 139 140C Initialization section 141 IF (FIRST_TIME .OR. IA.LE.0 .OR. range.GT.RGLAST) THEN 142C Find if cell or range have changed 143 samcel = .true. 144 if (first_time) then 145 samcel = .false. 146 else 147 do IX = 1,NX 148 do JX = 1,NX 149 IF (CELL(JX,IX) .NE. CELAST(JX,IX)) samcel = .false. 150 enddo 151 enddo 152 if (range .NE. RGLAST) samcel = .false. 153 endif 154 155C Cell initializations 156 if (.not.samcel) then 157C Store cell and range for comparison in subsequent calls 158 do IX = 1,NX 159 do JX = 1,NX 160 CELAST(JX,IX) = CELL(JX,IX) 161 enddo 162 enddo 163 RGLAST = range 164 FIRST_TIME = .false. 165 166C Notify to rangeR that CELL has changed 167 call mranger( 'CELL', CELL, range, NA, XA, NA, IAMOVE, 168 & IA, ISC, X0, NNA, MAXNNA ) 169 endif 170 171C Notify to rangeR that atoms have moved 172 call mranger( 'MOVE', CELL, range, NA, XA, NA, IAMOVE, 173 & IA, ISC, X0, NNA, MAXNNA ) 174 175 endif 176 177C Find neighbours of atom IA 178 if (IA .GT. 0) 179 & call mranger( 'FIND', CELL, range, NA, XA, NA, IAMOVE, 180 & IA, ISC, X0, NNA, MAXNNA ) 181 182C Find neighbours of point centered at x0 183C (not used unless MODE='FIND') 184C The point x0 is introduced as a variable of the module 185 if (IA .EQ. 0) then 186 call mranger( 'FIND', CELL, range, NA, XA, NA, IAMOVE, 187 & IA, ISC, X0, NNA, MAXNNA ) 188 endif 189 190 end subroutine mneighb 191 192 subroutine mranger( mode, cell, range, na, xa, 193 & namove, iamove, ia0, isc, x0, 194 & nna, maxnna ) 195 196C ******************************************************************** 197C Finds the neighbours of an atom in a cell with periodic boundary 198C conditions. Alternatively, it finds the atoms within a sphere 199C centered at an arbitrary point. It also allows to update the atomic 200C positions one at a time, what is useful in Montecarlo simulations. 201C Written by J.M.Soler. Nov'96. 202C *********** INPUT ************************************************** 203C CHARACTER*4 MODE : MODE='CELL' => Initialize or reshape cell 204C MODE='MOVE' => Move atom(s) 205C MODE='FIND' => Find neighbours 206C REAL*8 CELL(NX,NX) : Unit cell vectors CELL(IXYZ,IVECT) 207C REAL*8 range : Maximum distance of neighbours required 208C INTEGER NA : Number of atoms 209C REAL*8 XA(NX,NA) : Atomic positions in cartesian coordinates 210C INTEGER NAMOVE : Number of atoms to be moved 211C (not used unless MODE='MOVE') 212C INTEGER IAMOVE(NAMOVE) : Index(es) of atom(s) moved (not used 213C unless MODE='MOVE' and 0<NAMOVE<NA) 214C INTEGER IA0 : Atom whose neighbours are needed. 215C If IA0=0, point X0 is used as origin instead 216C (not used unless MODE='FIND') 217C INTEGER ISC : Single-counting switch (0=No, 1=Yes). 218C If ISC=1, only neighbours with JA.LE.IA0 219C are included in JAN 220C (not used unless MODE='FIND' and IA0.NE.0) 221C REAL*8 X0(NX) : Origin from which atoms are to be found, 222C in cartesian coordinates. 223C (not used unless MODE='FIND' and IA0=0) 224C INTEGER MAXNNA : Size of arrays JAN, XIJ and R2IJ 225C (not used unless MODE='FIND') 226C *********** OUTPUT ************************************************* 227C REAL*8 CELL(NX,NX) : Unit cell vectors CELL(IXYZ,IVECT) 228C The output cell is generated only when input 229C CELL has zero volume 230C INTEGER NNA : Number of 'neighbour' atoms within range of 231C atom IA0 or position X0 (only for MODE='FIND') 232C INTEGER JAN(NNAin) : Atom-index of neighbours (only for MODE='FIND') 233C REAL*8 XIJ(NX,NNAin): Vectors from atom IA0 or point X0 to neighbours 234C in cartesian coordinates (only for MODE='FIND') 235C REAL*8 R2IJ(NNAin) : Squared distances to neighbours 236C (only for MODE='FIND') 237C *********** UNITS ************************************************** 238C Units of CELL, range, XA and X0 are arbitrary but must be equal. 239C All vectors in cartesian coordinates. 240C *********** SUBROUTINES USED *************************************** 241C DISMIN, RECCEL, VOLCEL 242C *********** BEHAVIOUR ********************************************** 243C This is a 'remembering' routine, that saves a single copy of required 244C information on the system. Therefore, it cannot be used 245C simultaneously for different cells or sets of atoms. 246C A call with MODE='CELL' is required if the cell is changed. This also 247C updates all atomic positions with no need of a MODE='MOVE' call. 248C A call with MODE='MOVE' is required if any atoms are moved without 249C reshaping the cell, before any subsequent 'FIND' calls. 250C The routine knows if it has not been ever called, so that initial 251C calls with MODE='CELL' and MODE='MOVE' are implicit but not required 252C If MODE='MOVE' and NAMOVE=NA, all the atomic positions are 253C reinitialized, and the list IAMOVE is not used. This is also 254C true if MODE='CELL', irrespective of the value of NAMOVE 255C This routine works always with periodic boundary conditions. 256C If periodic boundary conditions are not desired, you must either 257C define a CELL large enough to contain all the atoms without 258C 'interaction' (i.e. distances smaller than range) between 259C different cells, or make CELL vectors identical to zero, in which 260C case an appropiate CELL is generated automatically by rangeR. 261C The size of this cell is determined by parameters DXMARG and DXRANG 262C below, and should be safe for small atomic displacements, but notice 263C that, if atoms move too much after the cell is generated (i.e. after 264C the last MODE='CELL' call), spureous 'interactions' may occur. 265C Do not make CELL extremely large to avoid this, since internal 266C array memory increases with cell volume. 267C If the number of neighbour atoms found is larger than the size of 268C the arrays JAN, XIJ and R2IJ, i.e. if NNAout > NNAin, these arrays 269C are filled only up to their size NNAin. With dynamic memory 270C allocation, this allows to find first the required array sizes 271C and then find the neighbours. Notice however that no warning is 272C given, so that you should always check that NNAout.LE.NNAin. 273C Different ranges can be used for different atoms or origins, but for 274C good performance, the largest range should be used in the initial 275C call (with MODE='CELL'). 276C There are no limitations regarding cell shape or size. The range may 277C be larger than the cell size, in which case many 'images' of the 278C same atom will be included in the neighbour list, with different 279C interatomic vectors XIJ and distances R2IJ. 280C The atom IA0 itself is included in the neighbour list, with zero 281C distance. You have to discard it if you want so. 282C CPU time and memory scale linearly with the number of atoms, for 283C sufficiently large numbers. 284C This is not an extremely optimized routine. The emphasis has been 285C rather put on functionality. It is not intended to vectorize or 286C parallelize well either. However, it is expected to be very 287C reasonably efficient on scalar machines, since most of the time 288C will be spent on the relatively simple and optimized 'search 289C section' at the end. 290C Code from beginning of loop 120 to label 130 was designed to exclude 291C neighbour cells which are inside a parallelepiped containing the 292C range sphere, but outside this sphere. This takes a nonneglegible 293C CPU time compared to one 'FIND' call for each atom (i.e. a MD step), 294C but may be worth if CELL is fixed or changes rarely compared to the 295C atomic displacements. But this is only true when CELL is strongly 296C nonorthorhombic (like fcc, bcc or hex) or when parameter NCR.GE.3 297C If you have doubts, don't worry: all this is normally irrelevant 298C unless you do Parrinello-Rahman dynamics. 299C *********** USAGE ************************************************** 300C Example of usage for a molecular dynamics simulation: 301C DIMENSION CELL(3,3),JAN(MAXNNA),R2IJ(MAXNNA),XIJ(3,MAXNNA),XA(3,NA) 302C Define CELL and initial positions XA 303C DO ITER = 1,NITER (Molecular dynamics iteration) 304C CALL rangeR('MOVE',3,CELL,range,NA,XA,NA,0,0,1,X0, 305C NNA,JAN,XIJ,R2IJ) 306C Initialize to zero all atomic forces FA(IX,IA) 307C DO IA = 1,NA (Loop on atoms) 308C NNA = MAXNNA 309C CALL rangeR('FIND',3,CELL,range,NA,XA,0,0,IA,1,X0, 310C NNA,JAN,XIJ,R2IJ) 311C IF (NNA.GT.MAXNNA) STOP 'Parameter MAXNNA too small' 312C DO IN = 1,NNA (Loop on neighbours of IA) 313C JA = JAN(IN) (Atomic index of neighbour) 314C RIJ = SQRT(R2IJ(IN)) (Interatomic distance) 315C IF (RIJ .GT. 1.D-12) THEN (Discard atom itself) 316C Find interatomic force FIJ( RIJ ) 317C DO IX = 1,3 318C FA(IX,IA) = FA(IX,IA) - FIJ * XIJ(IX,IN) / RIJ 319C FA(IX,JA) = FA(IX,JA) + FIJ * XIJ(IX,IN) / RIJ 320C ENDDO 321C ENDIF 322C ENDDO 323C ENDDO 324C Move atomic positions XA (Molecular dynamics step) 325C ENDDO 326C 327C Example of usage for a Montecarlo simulation: 328C DIMENSION CELL(3,3),JAN(MAXNNA),R2IJ(MAXNNA),XIJ(3,MAXNNA), 329C XA(3,NA),XNEW(3) 330C Define CELL and initial positions XA 331C CALL rangeR('CELL',3,CELL,range,NA,XA,0,0,0,0,XNEW, 332C NNA,JAN,XIJ,R2IJ) 333C Find here each atom's interaction energy EA(IA) 334C DO ITER = 1,NITER (Montecarlo iteration) 335C Choose atom moved IA and its trial new position XNEW 336C NNA = MAXNNA 337C CALL rangeR('FIND',3,CELL,range,NA,XA,0,0,0,0,XNEW, 338C NNA,JAN,XIJ,R2IJ) 339C IF (NNA.GT.MAXNNA) STOP 'Parameter MAXNNA too small' 340C EANEW = 0.D0 341C DO IN = 1,NNA (Loop on neighbours of IA) 342C JA = JAN(IN) (Atomic index of neighbour) 343C IF (JA.NE.IA) (Discard atom itself) 344C RIJ = SQRT(R2IJ(IN)) (Interatomic distance) 345C EANEW = EANEW + VIJ(RIJ) (Add interaction energy with JA) 346C ENDIF 347C ENDDO 348C IF (EXP(-(EANEW-EA(IA))/TEMP).GT.RAND()) THEN (Move atom) 349C EA(IA) = EANEW 350C DO IX = 1,3 351C XA(IX,IA) = XNEW(IX) 352C ENDDO 353C CALL rangeR('MOVE',3,CELL,range,NA,XA,1,IA,0,0,XNEW, 354C NNA,JAN,XIJ,R2IJ) 355C ENDIF 356C ENDDO 357C *********** ALGORITHM ********************************************** 358C The unit cell is divided in 'mesh-cells', and a list of atoms in 359C each cell is stored. To look for the neighbours of an atom, the 360C distances to the atoms in the neighbour cells are calculated and 361C compared to the range. 362C The list of atoms within a cell is stored as an 'ordered list', of 363C the kind so popular in C, using pointers from an atom to the next. 364C To deal with periodic boundary conditions, the mesh is 'extended' 365C on each side of the unit cell, and an index from the extended to 366C the unextended meshes is stored. In the extended mesh, the 367C 'index-shifts' and the vector distances from a mesh point (within 368C the unit cell) to its neighbour mesh points are independent of 369C the point, and they can thus be stored only once. 370C To calculate the vector from an atom to its neighbours, the 371C position of each atom relative to its mesh cell is also stored. 372C Thus, the vector between two atoms is the vector between their 373C cells plus the difference between their positions relative to 374C their cells. 375C *********** LANGUAGE *********************************************** 376C Illegal ANSI Fortran77 features used: IMPLICIT NONE, INCLUDE 377C Illegal ANSI Fortran90 features used: none. 378C *********** HISTORY AND CHANGES ************************************ 379C This is an improved version of routine NEIGHB, first written in C by 380C J.M.Soler in March 1995 and later translated to Fortran in Nov'96 381C The added capabilities of rangeR over NEIGHB are: 382C - Using a variable space dimension 383C - Finding the atoms within a sphere centered at an arbitrary point 384C - Updating the position of a single atom 385C - Finding the number of neighbours without any arrays 386C - Automatic cell generation for nonperiodic boundary conditions 387C All the algorithms and coding developed by J.M.Soler. November 1996. 388C ******************************************************************** 389 390C 391C Modules 392C 393 use precision 394 use alloc 395 396 IMPLICIT NONE 397 398C Argument types and dimensions 399 CHARACTER MODE*4 400 INTEGER NA, NAMOVE, NNA, MAXNNA 401 INTEGER IA0, IAMOVE(*), ISC 402 real(dp) :: 403 & CELL(*), range, X0(NX), XA(NX,NA) 404 405C NCR is the ratio between range radius and mesh-planes distance. 406C It fixes the size (and number) of mesh cells. 407C Recommended values are between 1 and 3 408 INTEGER NCR 409 PARAMETER ( NCR = 2 ) 410 411C DXMARG and DXRANG are used for automatic CELL generation 412C DXMARG is the minimum margin relative to coordinate range 413C DXRANG is the minimum margin relative to range 414C EPS is a small number to be subtracted from 1 415 real(dp) :: DXMARG, DXRANG, EPS 416 PARAMETER ( DXMARG = 0.1D0 ) 417 PARAMETER ( DXRANG = 1.0D0 ) 418 PARAMETER ( EPS = 1.D-14 ) 419 420C Variable-naming hints: 421C Character(s) indicates 422C A Atom 423C D Difference 424C E Extended (mesh) 425C I,J Index 426C M Mesh cell or point (lower vertex of mesh cell) 427C N Neighbour or Number (if first character) 428C R distance or vector modulus 429C V Vertex 430C X coordinate or basis vector 431C 1 lower bound 432C 2 upper bound or square 433 434C Internal functions, variables and arrays 435C REAL*8 CELMSH(MX*MX) Mesh-cell vectors 436C REAL*8 DMX(MX) In-cell atomic position in mesh coordinates 437C REAL*8 DPLANE Distance between lattice or mesh planes 438C REAL*8 DRM Minimum distance between two mesh cells 439C REAL*8 DX(3) Vector between two atoms 440C REAL*8 DX0M(MX) Origin position within mesh cell 441C REAL*8 DXAM(MX) Atom position within mesh cell 442C REAL*8 DXM(MX) Minimum vector between two mesh cells 443C REAL*8 DXMARG Parameter defined above 444C REAL*8 DXNM(MX,MAXNM) Cartesian vector between neighbour mesh points 445C REAL*8 DXRANG Parameter defined above 446C REAL*8 EPS Parameter defined above 447C INTEGER IA Atom index 448C INTEGER IAM Atom-to-move index 449C INTEGER IA1M(NM) Pointer to first atom in mesh cell 450C INTEGER IANEXT(NA) Pointer to next atom in mesh cell 451C INTEGER IAPREV(NA) Pointer to previous atom in mesh cell 452C INTEGER IDNM(MAXNM) Index-distance between neighbour mesh points 453C INTEGER IM Mesh index 454C INTEGER IEM Extended-mesh index 455C INTEGER IEMA(NA) Extended-mesh index of atoms 456C INTEGER I1EMX(MX) Minimum vaue of extended-mesh-coordinate indices 457C INTEGER I2EMX(MX) Maximum vaue of extended-mesh-coordinate indices 458C INTEGER I1MX(MX) Minimum vaue of mesh-coordinate indices 459C INTEGER I2MX(MX) Maximum vaue of mesh-coordinate indices 460C INTEGER IMX(MX) Mesh-cell index for each mesh vector 461C INTEGER IN Neighbour-mesh-cell index 462C INTEGER I1NX(MX) Minimum neighbour-cell-coordinate indices 463C INTEGER I2NX(MX) Maximum neighbour-cell-coordinate indices 464C INTEGER INX(MX) Neighbour-cell-coordinate indices 465C LOGICAL INSIDE Are two mesh cells within each other's range? 466C INTEGER IMESH(NEM) Correspondence between extended and normal mesh 467C INTEGER IV Vertex index 468C INTEGER IVX Vertex coordinate index 469C INTEGER IX Cartesian coordinte index 470C INTEGER IXX Double cartesian coordinte index 471C INTEGER JA Atom index 472C INTEGER JEM Extended-mesh index 473C INTEGER J1NX(MX) Minimum vertex-coordinate indices 474C INTEGER J2NX(MX) Maximum vertex-coordinate indices 475C INTEGER JX Cartesian coordinte index 476C LOGICAL MOVALL Move all atoms? 477C INTEGER NAM Number of atoms to move 478C INTEGER NCR Parameter defined above 479C INTEGER NEM Number of extended-mesh cells 480C INTEGER NEMX(MX) Extended-mesh cells in each mesh direction 481C INTEGER NM Number of mesh cells 482C INTEGER NMX(MX) Mesh cells in each mesh direction 483C INTEGER NNM Number of neighbour mesh cells 484C INTEGER NNMMAX Maximum number of neighbour mesh cells 485C INTEGER NNX(MX) Neighbour-cell ranges 486C LOGICAL NULCEL Null cell? 487C REAL*8 R2 Squared distance between two atoms 488C REAL*8 range2 Square of range 489C REAL*8 RCELL(MX*MX) Reciprocal cell vectors 490C RECCEL() Finds reciprocal lattice vectors 491C REAL*8 RNGMAX Maximum range 492C REAL*8 RMCELL(MX*MX) Reciprocal mesh-cell vectors 493C REAL*8 Rrange Slightly reduced range 494C REAL*8 XDIFF Range of atom coordinates 495C REAL*8 XMAX Maximum atom coordinate 496C REAL*8 XMIN Minimum atom coordinate 497 498 INTEGER 499 & IA, IAM, IEM, IM, 500 & IN, IX, IXX, JA, JEM, JM, JX, 501 & NAM, NM, NEM, NNM, NNMMAX 502c 503c Auxiliary variable to avoid compiler warnings 504c 505 integer j_aux 506 507 real(dp) 508 & DISMIN, DDOT, DPLANE, 509 & R2, range2, RNGMAX, Rrange, 510 & XDIFF, XMARG, XMAX, XMIN 511 512 513 logical 514 & INSIDE, MOVALL, NULCEL 515 516 external DISMIN, DDOT 517 518 save 519 & IAM, IEM, IM, 520 & NEM, NM, NNM, range2, RNGMAX, Rrange 521 522C Allocate local memory - check for change in number of atoms 523C and if there has been one then re-initialise 524 if (NA.gt.MAXNA) then 525 if (MAXNA.eq.-1) nullify(IANEXT,IAPREV,IEMA,DXAM) 526 call re_alloc( IANEXT, 1, NA, 'ianext', 'neighbour' ) 527 call re_alloc( IAPREV, 1, NA, 'iaprev', 'neighbour' ) 528 call re_alloc( IEMA, 1, NA, 'iema', 'neighbour' ) 529 call re_alloc( DXAM, 1, NX, 1, NA, 'dxam', 'neighbour' ) 530 MAXNA = NA 531 endif 532 533C Cell-mesh initialization section 534 IF (MODE.EQ.'CELL' .OR. MODE.EQ.'cell' .OR. 535 & range.GT.RNGMAX) THEN 536 537C Start time counter (this is for debugging) 538* CALL TIMER( 'rangeR1', 1 ) 539 540C Store range for comparison in subsequent calls 541 RNGMAX = range 542 543C Reduce the range slitghtly to avoid numerical-roundoff 544C ambiguities 545 Rrange = range * (1.D0 - EPS) 546 range2 = Rrange**2 547 548C Check if CELL must be generated automatically 549 NULCEL = .true. 550 DO 20 IXX = 1,NX*NX 551 IF (CELL(IXX) .NE. 0.D0) NULCEL = .false. 552 20 CONTINUE 553 IF (NULCEL) THEN 554 DO 40 IX = 1,NX 555C Find atom position bounds 556 XMIN = 1.D30 557 XMAX = -1.D30 558 DO 30 IA = 1,NA 559 XMIN = MIN( XMIN, XA(IX,IA) ) 560 XMAX = MAX( XMAX, XA(IX,IA) ) 561 30 CONTINUE 562C Determine 'cell margins' to prevent intercell interactions 563 XDIFF = XMAX - XMIN 564 XMARG = MAX( range*DXRANG, XDIFF*DXMARG ) 565C Define orthorrombic cell 566 IXX = IX + NX * (IX-1) 567 CELL(IXX) = XDIFF + 2.D0 * XMARG 568 40 CONTINUE 569 ENDIF 570 571C Find reciprocal cell vectors (not multiplied by 2*pi) 572 CALL RECCEL( NX, CELL, RCELL, 0 ) 573 574C Find number of mesh divisions 575 NM = 1 576 DO 50 IX = 1,NX 577 IXX = 1 + NX * (IX-1) 578 DPLANE = 1.D0 / SQRT(DDOT(NX,RCELL(IXX),1,RCELL(IXX),1) ) 579 NMX(IX) = 0.999D0 * DPLANE / (Rrange / NCR) 580 IF (NMX(IX) .LE. 0) NMX(IX) = 1 581 NM = NM * NMX(IX) 582 50 CONTINUE 583 584C Find mesh-cell vectors 585 IXX = 0 586 DO 70 IX = 1,NX 587 DO 60 JX = 1,NX 588 IXX = IXX + 1 589 CELMSH(IXX) = CELL(IXX) / NMX(IX) 590 RMCELL(IXX) = RCELL(IXX) * NMX(IX) 591 60 CONTINUE 592 70 CONTINUE 593 594C Find index-range of neighbour mesh cells and of extended mesh 595 NNM = 1 596 NEM = 1 597 DO 80 IX = 1,NX 598 IXX = 1 + NX * (IX-1) 599 DPLANE = 1.D0 / SQRT(DDOT(NX,RMCELL(IXX),1,RMCELL(IXX),1) ) 600 NNX(IX) = Rrange / DPLANE + 1 601 J1NX(IX) = 0 602 J2NX(IX) = 1 603 I1NX(IX) = - NNX(IX) 604 I2NX(IX) = + NNX(IX) 605 I1MX(IX) = 0 606 I2MX(IX) = NMX(IX) - 1 607 I1EMX(IX) = - NNX(IX) 608 I2EMX(IX) = NMX(IX) + NNX(IX) - 1 609 NEMX(IX) = NMX(IX) + 2 * NNX(IX) 610 NNM = NNM * (1+2*NNX(IX)) 611 NEM = NEM * NEMX(IX) 612 80 CONTINUE 613 614C Allocate arrays whose dimensions are now known 615 if (NM.gt.MAXNM) then 616 if (MAXNM.eq.-1) nullify(IA1M) 617 call re_alloc( IA1M, 1, NM, 'ia1m', 'neighbour' ) 618 MAXNM = NM 619 endif 620 if (NNM.gt.MAXNNM) then 621 if (MAXNNM.eq.-1) nullify(IDNM) 622 call re_alloc( IDNM, 1, NNM, 'idnm', 'neighbour' ) 623 call re_alloc( DXNM, 1, NX, 1, NNM, 'dxnm', 'neighbour' ) 624 MAXNNM = NNM 625 endif 626 627 if (NEM.gt.MAXNEM) then 628 if (MAXNEM.eq.-1) nullify(IMESH) 629 call re_alloc( IMESH, 1, NEM, 'imesh', 'neighbour' ) 630 MAXNEM = NEM 631 endif 632 633C Find which mesh cells are actually within range 634 NNMMAX = NNM 635 NNM = 0 636 DO 170 IN = 1,NNMMAX 637 j_aux = in 638 CALL INDARR( -1, NX, I1NX, I2NX, INX, 1, j_aux ) 639 INSIDE = .true. 640C From here to label 130 is generally not worth unless CELL is 641C very nonorthorrombic (like fcc, bcc or hex) and changes rarely 642* DO 120 IV = 1,2**NX 643* j_aux = iv 644* CALL INDARR( -1, NX, J1NX, J2NX, IVX, 1, j_aux ) 645* DO 100 IX = 1,NX 646* DXM(IX) = 0.D0 647* DO 90 JX = 1,NX 648* IXX = IX + NX * (JX-1) 649* DXM(IX) = DXM(IX) + CELMSH(IXX) * (INX(JX)+IVX(JX)) 650* 90 CONTINUE 651* 100 CONTINUE 652* DRM = DISMIN( NX, CELMSH, DXM ) 653* IF (DRM .LT. Rrange) THEN 654* INSIDE = .true. 655* GOTO 130 656* ENDIF 657* 120 CONTINUE 658* INSIDE = .false. 659* 130 CONTINUE 660 IF (INSIDE) THEN 661 NNM = NNM + 1 662C IDNM is the extended-mesh-index distance between 663C neighbour mesh cells 664 IDNM(NNM) = INX(NX) 665 DO 140 IX = NX-1,1,-1 666 IDNM(NNM) = INX(IX) + NEMX(IX) * IDNM(NNM) 667 140 CONTINUE 668C DXNM is the vector distance between neighbour mesh cells 669 DO 160 IX = 1,NX 670 DXNM(IX,NNM) = 0.D0 671 DO 150 JX = 1,NX 672 IXX = IX + NX * (JX-1) 673 DXNM(IX,NNM) = DXNM(IX,NNM) + CELMSH(IXX) * INX(JX) 674 150 CONTINUE 675 160 CONTINUE 676 ENDIF 677 170 CONTINUE 678 679C Find correspondence between extended and reduced (normal) meshes 680 do IEM = 1,NEM 681 j_aux = iem 682 CALL INDARR( -1, NX, I1EMX, I2EMX, IMX, 1, j_aux ) 683 CALL INDARR( +1, NX, I1MX, I2MX, IMX, 1, IM ) 684 IMESH(IEM) = IM 685 enddo 686 687C Stop time counter 688* CALL TIMER( 'rangeR1', 2 ) 689 690C Set 'move all atoms' switch 691 MOVALL = .true. 692 ELSE 693 MOVALL = .false. 694 ENDIF 695C End of cell initialization section 696 697C Atom-positions (relative to mesh) initialization section 698 IF (MODE.EQ.'MOVE' .OR. MODE.EQ.'move' .OR. MOVALL) THEN 699 IF (NAMOVE .EQ. NA) MOVALL = .true. 700 701C Start time counter 702* CALL TIMER( 'rangeR2', 1 ) 703 704 IF (MOVALL) THEN 705 NAM = NA 706C Initialize 'atoms in mesh-cell' lists 707 do IA = 1,NA 708 IANEXT(IA) = 0 709 IAPREV(IA) = 0 710 enddo 711 do IM = 1,NM 712 IA1M(IM) = 0 713 enddo 714 ELSE 715 NAM = NAMOVE 716 ENDIF 717 718C Loop on moved atoms 719 DO 240 IAM = 1,NAM 720 721C Select atom to move 722 IF (MOVALL) THEN 723 IA = IAM 724 ELSE 725 IA = IAMOVE(IAM) 726C Supress atom from its previous mesh-cell 727 JA = IAPREV(IA) 728 IF (JA.NE.0) IANEXT(JA) = IANEXT(IA) 729 JA = IANEXT(IA) 730 IF (JA.NE.0) IAPREV(JA) = IAPREV(IA) 731 IEM = IEMA(IA) 732 IM = IMESH(IEM) 733 IF (IA1M(IM) .EQ. IA) IA1M(IM) = JA 734 ENDIF 735 736C Find mesh-cell in which atom is 737 DO 220 IX = 1,NX 738 IXX = 1 + NX * (IX-1) 739 DMX(IX) = DDOT(NX,RMCELL(IXX),1,XA(1,IA),1) 740 IMX(IX) = INT( DMX(IX) + 1000.D0 ) - 1000 741 DMX(IX) = DMX(IX) - IMX(IX) 742 IMX(IX) = MOD( IMX(IX) + 1000 * NMX(IX), NMX(IX) ) 743 220 CONTINUE 744 CALL INDARR( +1, NX, I1EMX, I2EMX, IMX, 1, IEM ) 745 CALL INDARR( +1, NX, I1MX, I2MX, IMX, 1, IM ) 746 IEMA(IA) = IEM 747 748C Put atom first in its new mesh-cell 749 JA = IA1M(IM) 750 IF (JA .NE. 0) IAPREV(JA) = IA 751 IANEXT(IA) = JA 752 IA1M(IM) = IA 753 754C Find atomic position relative to mesh 755 DO 230 IX = 1,NX 756 DXAM(IX,IA) = 0.D0 757 DO 225 JX = 1,NX 758 IXX = IX + NX * (JX-1) 759 DXAM(IX,IA) = DXAM(IX,IA) + CELMSH(IXX) * DMX(JX) 760 225 CONTINUE 761 230 CONTINUE 762 763 240 CONTINUE 764 765C Stop time counter 766C CALL TIMER( 'rangeR2', 2 ) 767 ENDIF 768C End of atom-positions initialization section 769 770C Search section 771 IF (MODE.EQ.'FIND' .OR. MODE.EQ.'find') THEN 772 Rrange = range * (1.D0 - EPS) 773 range2 = range**2 774 775C Find the mesh cell of the center of the sphere 776 IF (IA0.LE.0) THEN 777C Find mesh cell of position X0 778 DO 250 IX = 1,NX 779 IXX = 1 + NX * (IX-1) 780 DMX(IX) = DDOT(NX,RMCELL(IXX),1,X0,1) 781 IMX(IX) = INT( DMX(IX) + 1000.D0 ) - 1000 782 DMX(IX) = DMX(IX) - IMX(IX) 783 IMX(IX) = MOD( IMX(IX) + 1000 * NMX(IX), NMX(IX) ) 784 250 CONTINUE 785 CALL INDARR( +1, NX, I1EMX, I2EMX, IMX, 1, IEM ) 786 DO 270 IX = 1,NX 787 DX0M(IX) = 0.D0 788 DO 260 JX = 1,NX 789 IXX = IX + NX * (JX-1) 790 DX0M(IX) = DX0M(IX) + CELMSH(IXX) * DMX(JX) 791 260 CONTINUE 792 270 CONTINUE 793 ELSE 794C Find mesh cell of atom IA0 795 IEM = IEMA(IA0) 796 DO 280 IX = 1,NX 797 DX0M(IX) = DXAM(IX,IA0) 798 280 CONTINUE 799 ENDIF 800 801C Loop on neighbour mesh cells and on the atoms within them 802C This is usually the only time-consuming loop 803 NNA = 0 804 do IN = 1,NNM 805 JEM = IEM + IDNM(IN) 806 JM = IMESH(JEM) 807C Loop on atoms of neighbour cell. 808C Try first atom in this mesh-cell 809 JA = IA1M(JM) 810 300 CONTINUE 811 IF (JA .NE. 0) THEN 812C Check that single-counting exclusion does not apply 813 IF (IA0.LE.0 .OR. ISC.EQ.0 .OR. JA.LE.IA0) THEN 814C Find vector and distance to atom JA 815 R2 = 0.0d0 816 do IX = 1,NX 817 DX(IX) = DXNM(IX,IN) + DXAM(IX,JA) - DX0M(IX) 818 R2 = R2 + DX(IX)**2 819 enddo 820C Check if atom JA is within range 821 IF (R2 .LE. range2) THEN 822 NNA = NNA + 1 823C Check that array arguments are not overflooded 824 IF (NNA .GT. MAXNNA) THEN 825 call sizeup_neighbour_arrays( MAXNNA+NA ) 826 ENDIF 827 JAN(NNA) = JA 828 do IX = 1,NX 829 XIJ(IX,NNA) = DX(IX) 830 enddo 831 R2IJ(NNA) = R2 832 ENDIF 833 ENDIF 834C Take next atom in this mesh-cell and go to begining of loop 835 JA = IANEXT(JA) 836 GOTO 300 837 ENDIF 838 enddo 839 ENDIF 840C End of search section 841 return 842 end subroutine mranger 843 844 subroutine reccel( N, A, B, IOPT ) 845 846C CALCULATES RECIPROCAL LATTICE VECTORS B. 847C THEIR PRODUCT WITH DIRECT LATTICE VECTORS A IS 1 (IF IOPT=0) OR 848C 2*PI (IF IOPT=1). N IS THE SPACE DIMENSION. 849C WRITTEN BY J.M.SOLER. 850 851 integer :: n, iopt 852 real(dp):: A(N,N),B(N,N), c, ci 853 854 integer :: i 855 856 C=1.D0 857 IF (IOPT.EQ.1) C=2.D0*ACOS(-1.D0) 858 859 IF (N .EQ. 1) THEN 860 B(1,1) = C / A(1,1) 861 ELSEIF (N .EQ. 2) THEN 862 C = C / (A(1,1)*A(2,2) - A(1,2)*A(2,1)) 863 B(1,1) = A(2,2)*C 864 B(1,2) = (-A(2,1))*C 865 B(2,1) = (-A(1,2))*C 866 B(2,2) = A(1,1)*C 867 ELSEIF (N .EQ. 3) THEN 868 B(1,1)=A(2,2)*A(3,3)-A(3,2)*A(2,3) 869 B(2,1)=A(3,2)*A(1,3)-A(1,2)*A(3,3) 870 B(3,1)=A(1,2)*A(2,3)-A(2,2)*A(1,3) 871 B(1,2)=A(2,3)*A(3,1)-A(3,3)*A(2,1) 872 B(2,2)=A(3,3)*A(1,1)-A(1,3)*A(3,1) 873 B(3,2)=A(1,3)*A(2,1)-A(2,3)*A(1,1) 874 B(1,3)=A(2,1)*A(3,2)-A(3,1)*A(2,2) 875 B(2,3)=A(3,1)*A(1,2)-A(1,1)*A(3,2) 876 B(3,3)=A(1,1)*A(2,2)-A(2,1)*A(1,2) 877 DO 20 I=1,3 878 CI=C/(A(1,I)*B(1,I)+A(2,I)*B(2,I)+A(3,I)*B(3,I)) 879 B(1,I)=B(1,I)*CI 880 B(2,I)=B(2,I)*CI 881 B(3,I)=B(3,I)*CI 882 20 CONTINUE 883 ELSE 884 call die('RECCEL: NOT PREPARED FOR N>3') 885 ENDIF 886 end subroutine reccel 887 888 subroutine indarr( IOPT, ND, I1, I2, I, J1, J ) 889 890C ******************************************************************** 891C Finds the global index in a multidimensional array from the idexes 892C in each dimension, or viceversa (the first is an explicit solution 893C of the standard index-resolution problem that the compiler solves 894C each time an array element is referenced). 895C Written by J.M.Soler. Nov'96. 896C *********** INPUT ************************************************** 897C INTEGER IOPT : IOPT>0 => from I to J. IOPT<0 => from J to I. 898C INTEGER ND : Number of array dimensions (indexes) 899C INTEGER I1(ND) : Minimum value of array indexes 900C INTEGER I2(ND) : Maximum value of array indexes (i.e. ARRAY(I1:I2)) 901C INTEGER J1 : Minimum value of global index (typically 1) 902C *********** INPUT OR OUTPUT (DEPENDING OF IOPT) ******************** 903C INTEGER I(ND) : Array indexes in each dimension 904C INTEGER J : Global index (first dimension increasing fastest) 905C *********** BEHAVIOUR ********************************************** 906C Indexes I() are taken as periodic, i.e. their modulus I2(ID)-I1(ID)+1 907C is taken before using them. This simplifies its use as indexes of a 908C mesh with periodic boundary conditions. This modulus operation is 909C also done with J, so that the output I() are always within range. 910C If IOPT=0, nothing is done. 911C *********** USAGE ************************************************** 912C Sample usage to find the Laplacian of a function defined in a mesh 913C with periodic boundary conditions in a space of variable dimension 914C SUBROUTINE LAPLACIAN( ND, N, DX, F, FLAPL ) 915C PARAMETER (MAXD = 3) 916C DIMENSION N(ND),DX(ND),F(*),FLAPL(*),I1(MAXD),I2(MAXD),I(MAXD) 917C NMESH = 1 918C DO ID = 1,ND 919C I1(ID) = 1 920C I2(ID) = N(ID) 921C NMESH = NMESH * N(ID) 922C ENDDO 923C DO IMESH = 1, NMESH 924C CALL INDARR( -1, ND, I1, I2, I, 1, IMESH ) 925C FLAPL(IMESH) = 0. 926C DO ID = 1,ND 927C DO K = -1,1,2 928C I(ID) = I(ID) + K 929C CALL INDARR( +1, ND, I1, I2, I, 1, JMESH ) 930C FLAPL(IMESH) = FLAPL(IMESH) + F(JMESH) / DX(ID)**2 931C I(ID) = I(ID) - K 932C ENDDO 933C FLAPL(IMESH) = FLAPL(IMESH) - 2. * F(IMESH) / DX(ID)**2 934C ENDDO 935C ENDDO 936C END 937C *********** LANGUAGE *********************************************** 938C Illegal ANSI Fortran77 features used: IMPLICIT NONE 939C Illegal ANSI Fortran90 features used: none. 940C ******************************************************************** 941 942C Next line is non-standard but may be supressed 943 IMPLICIT NONE 944 INTEGER ND 945 INTEGER I(ND), I1(ND), I2(ND), ID, IOPT, J, J1, K, N 946 947 IF (IOPT .GT. 0) THEN 948 J = 0 949 DO 10 ID = ND,1,-1 950 N = I2(ID) - I1(ID) + 1 951 K = I(ID) - I1(ID) 952 K = MOD( K + 1000 * N, N ) 953 J = K + N * J 954 10 CONTINUE 955 J = J + J1 956 ELSEIF (IOPT .LT. 0) THEN 957 K = J - J1 958 DO 20 ID = 1,ND 959 N = I2(ID) - I1(ID) + 1 960 I(ID) = I1(ID) + MOD( K, N ) 961 K = K / N 962 20 CONTINUE 963 ENDIF 964 965 end subroutine indarr 966! 967 subroutine sizeup_neighbour_arrays(n) 968 implicit none 969 integer, intent(in) :: n 970 971 ! Makes sure that the neighbour arrays are at 972 ! least of size n 973 if (.not. pointers_allocated) then 974 nullify(jan) 975 nullify(r2ij) 976 nullify(xij) 977! Dimension arrays to initial size n 978 maxnna = n 979 call re_alloc( jan, 1, maxnna,'jan', 'neighbour' ) 980 call re_alloc( r2ij, 1, maxnna, 'r2ij', 'neighbour' ) 981 call re_alloc( xij, 1, 3, 1, maxnna, 'xij', 'neighbour' ) 982 pointers_allocated = .true. 983 else 984 if (n > maxnna) then 985 maxnna = n 986 call re_alloc( jan, 1, maxnna,'jan', 'neighbour' ) 987 call re_alloc( r2ij, 1, maxnna, 'r2ij', 'neighbour' ) 988 call re_alloc( xij, 1, 3, 1, maxnna, 'xij', 'neighbour' ) 989 endif 990 endif 991 end subroutine sizeup_neighbour_arrays 992 993 subroutine reset_neighbour_arrays( ) 994 implicit none 995!!#ifdef DEBUG 996 call write_debug( ' PRE reset_neighbour_arrays' ) 997!!#endif 998 999 celast = 0.0_dp 1000 rglast = 0.0_dp 1001 1002 if (pointers_allocated) then 1003 call de_alloc( jan, 'jan', 'neighbour' ) 1004 call de_alloc( r2ij, 'r2ij', 'neighbour' ) 1005 call de_alloc( xij, 'xij', 'neighbour' ) 1006 pointers_allocated = .false. 1007 endif 1008 1009 if (maxna.gt.0) then 1010 call de_alloc( IANEXT, 'ianext', 'neighbour' ) 1011 call de_alloc( IAPREV, 'iaprev', 'neighbour' ) 1012 call de_alloc( IEMA, 'iema', 'neighbour' ) 1013 call de_alloc( DXAM, 'dxam', 'neighbour' ) 1014 maxna = -1 1015 endif 1016 1017 if (maxnm.gt.0) then 1018 call de_alloc( IA1M, 'ia1m', 'neighbour' ) 1019 maxnm = -1 1020 endif 1021 1022 if (maxnnm.gt.0) then 1023 call de_alloc( IDNM, 'idnm', 'neighbour' ) 1024 call de_alloc( DXNM, 'dxnm', 'neighbour' ) 1025 maxnnm = -1 1026 endif 1027 1028 if (maxnem.gt.0) then 1029 call de_alloc( IMESH, 'imesh', 'neighbour' ) 1030 maxnem = -1 1031 endif 1032 1033!!#ifdef DEBUG 1034 call write_debug( ' POS reset_neighbour_arrays' ) 1035!!#endif 1036 end subroutine reset_neighbour_arrays 1037 1038 end module neighbour 1039 1040