1* 2* $Id$ 3* 4 subroutine incell2(ni,r2,r1) 5 6* ======================================================= 7* This routine places the ions inside the cell defined 8* by the lattice vectors centered at zero. 9* ======================================================= 10 implicit none 11 integer ni 12 real*8 r2(3,*),r1(3,*) 13 14* **** Local variables defined **** 15 real*8 fa1,fa2,fa3 16 real*8 a(3,3),b(3,3),volume 17 integer i,j 18 19* **** external functions **** 20 real*8 lattice_unita 21 external lattice_unita 22 23* ***** Determine the unit lattice vectors and distances ****** 24 do j=1,3 25 do i=1,3 26 a(i,j) = lattice_unita(i,j) 27 end do 28 end do 29 30 b(1,1) = a(2,2)*a(3,3) - a(3,2)*a(2,3) 31 b(2,1) = a(3,2)*a(1,3) - a(1,2)*a(3,3) 32 b(3,1) = a(1,2)*a(2,3) - a(2,2)*a(1,3) 33 b(1,2) = a(2,3)*a(3,1) - a(3,3)*a(2,1) 34 b(2,2) = a(3,3)*a(1,1) - a(1,3)*a(3,1) 35 b(3,2) = a(1,3)*a(2,1) - a(2,3)*a(1,1) 36 b(1,3) = a(2,1)*a(3,2) - a(3,1)*a(2,2) 37 b(2,3) = a(3,1)*a(1,2) - a(1,1)*a(3,2) 38 b(3,3) = a(1,1)*a(2,2) - a(2,1)*a(1,2) 39 volume = a(1,1)*b(1,1) 40 > + a(2,1)*b(2,1) 41 > + a(3,1)*b(3,1) 42 43 volume = 1.0d0/volume 44 call dscal(9,volume,b,1) 45 46!$OMP DO 47 do i =1,ni 48 49* *** Break the Ion positions into the a1, a2, and a3 components *** 50 fa1 = b(1,1) * r2(1,i) 51 > + b(2,1) * r2(2,i) 52 > + b(3,1) * r2(3,i) 53 54 fa2 = b(1,2) * r2(1,i) 55 > + b(2,2) * r2(2,i) 56 > + b(3,2) * r2(3,i) 57 58 fa3 = b(1,3) * r2(1,i) 59 > + b(2,3) * r2(2,i) 60 > + b(3,3) * r2(3,i) 61 62 63* **** Change the a1, a2 and a3 components to **** 64* **** make the ion be in the cell **** 65 66 23 IF (fa1 .GT. (0.5d0)) THEN 67* WRITE (*,*) 'a1>', I, R2A1, DA(1)/2.0d0 68 r2(1,i) = r2(1,i) - lattice_unita(1,1) 69 r2(2,i) = r2(2,i) - lattice_unita(2,1) 70 r2(3,i) = r2(3,i) - lattice_unita(3,1) 71 72 r1(1,i) = r1(1,i) - lattice_unita(1,1) 73 r1(2,i) = r1(2,i) - lattice_unita(2,1) 74 r1(3,i) = r1(3,i) - lattice_unita(3,1) 75 76 77 fa1 = b(1,1) * r2(1,i) 78 > + b(2,1) * r2(2,i) 79 > + b(3,1) * r2(3,i) 80 GO TO 23 81 ENDIF 82 83 24 IF (fa1 .LE. (-0.5d0)) THEN 84* WRITE (*,*) 'a1<', I, R2A1, DA(1)/2.0d0 85 r2(1,i) = r2(1,i) + lattice_unita(1,1) 86 r2(2,i) = r2(2,i) + lattice_unita(2,1) 87 r2(3,i) = r2(3,i) + lattice_unita(3,1) 88 89 r1(1,i) = r1(1,i) + lattice_unita(1,1) 90 r1(2,i) = r1(2,i) + lattice_unita(2,1) 91 r1(3,i) = r1(3,i) + lattice_unita(3,1) 92 93 94 fa1 = b(1,1) * r2(1,i) 95 > + b(2,1) * r2(2,i) 96 > + b(3,1) * r2(3,i) 97 GO TO 24 98 ENDIF 99 100 25 IF (fa2 .GT. (0.5d0)) THEN 101* WRITE (*,*) 'a2>', I, R2A2, DA(2)/2.0d0 102 r2(1,i) = r2(1,i) - lattice_unita(1,2) 103 r2(2,i) = r2(2,i) - lattice_unita(2,2) 104 r2(3,i) = r2(3,i) - lattice_unita(3,2) 105 106 r1(1,i) = r1(1,i) - lattice_unita(1,2) 107 r1(2,i) = r1(2,i) - lattice_unita(2,2) 108 r1(3,i) = r1(3,i) - lattice_unita(3,2) 109 110 111 fa2 = b(1,2) * r2(1,i) 112 > + b(2,2) * r2(2,i) 113 > + b(3,2) * r2(3,i) 114 GO TO 25 115 ENDIF 116 117 26 IF (fa2 .LE. (-0.5d0)) THEN 118* WRITE (*,*) 'a2<', I, R2A2, DA(2)/2.0d0 119 r2(1,i) = r2(1,i) + lattice_unita(1,2) 120 r2(2,i) = r2(2,i) + lattice_unita(2,2) 121 r2(3,i) = r2(3,i) + lattice_unita(3,2) 122 123 r1(1,i) = r1(1,i) + lattice_unita(1,2) 124 r1(2,i) = r1(2,i) + lattice_unita(2,2) 125 r1(3,i) = r1(3,i) + lattice_unita(3,2) 126 127 128 fa2 = b(1,2) * r2(1,i) 129 > + b(2,2) * r2(2,i) 130 > + b(3,2) * r2(3,i) 131 GO TO 26 132 ENDIF 133 134 135 27 IF (fa3 .GT. (0.5d0)) THEN 136* WRITE (*,*) 'a3>', i, R2A3, DA(3)/2.0d0 137 r2(1,i) = r2(1,i) - lattice_unita(1,3) 138 r2(2,i) = r2(2,i) - lattice_unita(2,3) 139 r2(3,i) = r2(3,i) - lattice_unita(3,3) 140 141 r1(1,i) = r1(1,i) - lattice_unita(1,3) 142 r1(2,i) = r1(2,i) - lattice_unita(2,3) 143 r1(3,i) = r1(3,i) - lattice_unita(3,3) 144 145 fa3 = b(1,3) * r2(1,i) 146 > + b(2,3) * r2(2,i) 147 > + b(3,3) * r2(3,i) 148 GO TO 27 149 ENDIF 150 151 28 IF (fa3 .LE. (-0.5d0)) THEN 152* WRITE (*,*) 'a3<', I, R2A3, DA(3)/2.0d0 153 r2(1,i) = r2(1,i) + lattice_unita(1,3) 154 r2(2,i) = r2(2,i) + lattice_unita(2,3) 155 r2(3,i) = r2(3,i) + lattice_unita(3,3) 156 157 r1(1,i) = r1(1,i) + lattice_unita(1,3) 158 r1(2,i) = r1(2,i) + lattice_unita(2,3) 159 r1(3,i) = r1(3,i) + lattice_unita(3,3) 160 161 fa3 = b(1,3) * r2(1,i) 162 > + b(2,3) * r2(2,i) 163 > + b(3,3) * r2(3,i) 164 GO TO 28 165 ENDIF 166 167 end do 168!$OMP END DO 169 170 return 171 end 172