1! This contains the algorithms needed to calculate the distance between atoms. 2! 3! Copyright (C) 2013, Joshua More and Michele Ceriotti 4! 5! Permission is hereby granted, free of charge, to any person obtaining a copy 6! of this software and associated documentation files (the "Software"), to deal 7! in the Software without restriction, including without limitation the rights 8! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell 9! copies of the Software, and to permit persons to whom the Software is 10! furnished to do so, subject to the following conditions: 11 12! The above copyright notice and this permission notice shall be included in 13! all copies or substantial portions of the Software. 14 15! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR 16! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 17! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE 18! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER 19! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, 20! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN 21! THE SOFTWARE. 22! 23! Functions: 24! vector_separation: Calculates the vector separating two atoms. 25! separation: Calculates the square distance between two vectors. 26! nearest_neighbours: Generates arrays to calculate the pairs of atoms within 27! a certain radius of each other. 28 29 MODULE DISTANCE 30 IMPLICIT NONE 31 32 CONTAINS 33 34 SUBROUTINE vector_separation(cell_h, cell_ih, ri, rj, rij, r2) 35 ! Calculates the vector separating two atoms. 36 ! 37 ! Note that minimum image convention is used, so only the image of 38 ! atom j that is the shortest distance from atom i is considered. 39 ! 40 ! Also note that while this may not work if the simulation 41 ! box is highly skewed from orthorhombic, as 42 ! in this case it is possible to return a distance less than the 43 ! nearest neighbour distance. However, this will not be of 44 ! importance unless the cut-off radius is more than half the 45 ! width of the shortest face-face distance of the simulation box, 46 ! which should never be the case. 47 ! 48 ! Args: 49 ! cell_h: The simulation box cell vector matrix. 50 ! cell_ih: The inverse of the simulation box cell vector matrix. 51 ! ri: The position vector of atom i. 52 ! rj: The position vector of atom j 53 ! rij: The vector separating atoms i and j. 54 ! r2: The square of the distance between atoms i and j. 55 56 DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN) :: cell_h 57 DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN) :: cell_ih 58 DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: ri 59 DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: rj 60 DOUBLE PRECISION, DIMENSION(3), INTENT(OUT) :: rij 61 DOUBLE PRECISION, INTENT(OUT) :: r2 62 63 INTEGER k 64 DOUBLE PRECISION, DIMENSION(3) :: sij 65 ! The separation in a basis where the simulation box 66 ! is a unit cube. 67 68 sij = matmul(cell_ih, ri - rj) 69 DO k = 1, 3 70 ! Finds the smallest separation of all the images of atom i and j 71 sij(k) = sij(k) - dnint(sij(k)) 72 ENDDO 73 rij = matmul(cell_h, sij) 74 r2 = dot_product(rij,rij) 75 76 END SUBROUTINE 77 78 SUBROUTINE separation(cell_h, cell_ih, ri, rj, r2) 79 ! Calculates the squared distance between two position vectors. 80 ! 81 ! Note that minimum image convention is used, so only the image of 82 ! atom j that is the shortest distance from atom i is considered. 83 ! 84 ! Also note that while this may not work if the simulation 85 ! box is highly skewed from orthorhombic, as 86 ! in this case it is possible to return a distance less than the 87 ! nearest neighbour distance. However, this will not be of 88 ! importance unless the cut-off radius is more than half the 89 ! width of the shortest face-face distance of the simulation box, 90 ! which should never be the case. 91 ! 92 ! Args: 93 ! cell_h: The simulation box cell vector matrix. 94 ! cell_ih: The inverse of the simulation box cell vector matrix. 95 ! ri: The position vector of atom i. 96 ! rj: The position vector of atom j 97 ! r2: The square of the distance between atoms i and j. 98 99 DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN) :: cell_h 100 DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN) :: cell_ih 101 DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: ri 102 DOUBLE PRECISION, DIMENSION(3), INTENT(IN) :: rj 103 DOUBLE PRECISION, INTENT(OUT) :: r2 104 105 INTEGER k 106 ! The separation in a basis where the simulation box 107 ! is a unit cube. 108 DOUBLE PRECISION, DIMENSION(3) :: sij 109 DOUBLE PRECISION, DIMENSION(3) :: rij 110 111 sij = matmul(cell_ih, ri - rj) 112 DO k = 1, 3 113 ! Finds the smallest separation of all the images of atom i and j 114 sij(k) = sij(k) - dnint(sij(k)) 115 ENDDO 116 rij = matmul(cell_h, sij) 117 r2 = dot_product(rij, rij) 118 119 END SUBROUTINE 120 121 SUBROUTINE nearest_neighbours(rn, natoms, atoms, cell_h, cell_ih, index_list, n_list) 122 ! Creates a list of all the pairs of atoms that are closer together 123 ! than a certain distance. 124 ! 125 ! This takes all the positions, and calculates which ones are 126 ! shorter than the distance rn. This creates two vectors, index_list 127 ! and n_list. index_list(i) gives the last index of n_list that 128 ! corresponds to a neighbour of atom i. 129 ! 130 ! 131 ! Args: 132 ! rn: The nearest neighbour list cut-off parameter. This should 133 ! be larger than the potential cut-off radius. 134 ! natoms: The number of atoms in the system. 135 ! atoms: A vector holding all the atom positions. 136 ! cell_h: The simulation box cell vector matrix. 137 ! cell_ih: The inverse of the simulation box cell vector matrix. 138 ! index_list: A array giving the last index of n_list that 139 ! gives the neighbours of a given atom. Essentially keeps 140 ! track of how many atoms neighbour a given atom. 141 ! n_list: An array giving the indices of the atoms that neighbour 142 ! the atom determined by index_list. Essentially keeps track 143 ! of which atoms neighbour a given atom. 144 145 DOUBLE PRECISION, INTENT(IN) :: rn 146 INTEGER, INTENT(IN) :: natoms 147 DOUBLE PRECISION, DIMENSION(:,:), INTENT(IN) :: atoms 148 DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN) :: cell_h 149 DOUBLE PRECISION, DIMENSION(3,3), INTENT(IN) :: cell_ih 150 INTEGER, DIMENSION(natoms), INTENT(OUT) :: index_list 151 INTEGER, DIMENSION(natoms*(natoms-1)/2), INTENT(OUT) :: n_list 152 153 INTEGER :: i, j 154 DOUBLE PRECISION r2 155 156 index_list(1) = 0 157 158 DO i = 1, natoms - 1 159 DO j = i + 1, natoms 160 CALL separation(cell_h, cell_ih, atoms(i,:), atoms(j,:), r2) 161 IF (r2 < rn*rn) THEN 162 ! We have found an atom that neighbours atom i, so the 163 ! i-th index of index_list is incremented by one, and a new 164 ! entry is added to n_list. 165 index_list(i) = index_list(i) + 1 166 n_list(index_list(i)) = j 167 ENDIF 168 ENDDO 169 index_list(i+1) = index_list(i) 170 ENDDO 171 172 END SUBROUTINE 173 174 END MODULE 175