1C Package tranal v.5.0 2C Calculations of radial distribution functions 3 program RDF 4 include "tranal.h" 5 character*64 FOUTRDF 6 parameter (MAXRDF=64,MAXGR=1024,NAAM=1024) 7 integer*8 IREP(MAXRDF),IRD(MAXRDF,MAXGR),JRD(MAXRDF,MAXGR), 8 +IRDF(MAXRDF,NAAM),IRDFI(MAXRDF,NAAM) 9 logical LCOMRDF1,LCOMRDF2 10 dimension TMP(MAXRDF) 11 namelist /RDFIN/ FOUTRDF,RDFCUT,NA,NAI,RMI,RMAX,NRDF,ME, 12 + LCOMRDF1,LCOMRDF2 13* defaults 14 LCOMRDF1=.false. 15 LCOMRDF2=.false. 16 RDFCUT=10. 17 NA=200 18 NAI=100 19 RMI=0 20 RMAX=5. 21 ME=-1 22* read and setup trajectory parameters 23 call SETTRAJ 24* read RDF computation parameters 25 read(*,RDFIN) 26 if(NRDF.gt.MAXRDF)stop 'Increase MAXRDF!!!' 27 if(NA.gt.NAAM.or.NAI.gt.NAAM)stop 'Increase NAAM!' 28 DRI=RMAX-RMI 29* initialization 30 BOXS=0. ! to accumulate average box sizes 31 BOYS=0. 32 BOZS=0. 33* read RDFs 34 N=0 35 if(LCOMRDF1)write(*,*) 36 +' For the first sites in the list, molecular COM is used' 37 if(LCOMRDF2)write(*,*) 38 +' For the second sites in the list, molecular COM is used' 39 write(*,*)' Trying to read ',NRDF,' groups for RDF calculations' 40 DO K = 1,NRDF 41 IE=-1 42 STR=TAKESTR(5,IE) 43 if(str(1:1).eq.'&')then 44 READ(str(2:24),*)IRP 45 if(IRP.gt.MAXGR) stop 'Increase MAXGR!!!' 46 IREP(K)=IRP 47 do IR=1,IRP 48 STR=TAKESTR(5,IE) 49 read(STR,*,err=99)I,J 50 IRD(K,IR)=I 51 JRD(K,IR)=J 52 end do 53 else 54 IREP(K)=1 55 read(STR,*,err=99)I,J 56 IRD(K,1)=I 57 JRD(K,1)=J 58 end if 59 write(*,*)' RDF No ',K,' equiv ',IREP(K) 60 if(LCOMRDF1.and.LCOMRDF2)then 61 do I=1,IREP(K) 62 write(*,'(2I6,2(3x,a6))') 63 + IRD(K,I),JRD(K,I),NAME(ITS(IRD(K,I))),NAME(ITS(JRD(K,I))) 64 end do 65 else if(LCOMRDF1.and..not.LCOMRDF2)then 66 do I=1,IREP(K) 67 write(*,'(2I6,2(3x,a6))') 68 + IRD(K,I),JRD(K,I),NAME(ITS(IRD(K,I))),NM(JRD(K,I)) 69 end do 70 else if(.not.LCOMRDF1.and.LCOMRDF2)then 71 do I=1,IREP(K) 72 write(*,'(2I6,2(3x,a6))') 73 + IRD(K,I),JRD(K,I),NM(IRD(K,I)),NAME(ITS(JRD(K,I))) 74 end do 75 else 76 do I=1,IREP(K) 77 write(*,'(2I6,2(3x,a4))') 78 + IRD(K,I),JRD(K,I),NM(IRD(K,I)),NM(JRD(K,I)) 79 end do 80 end if 81 END DO 82 go to 101 83 99 STR=TAKESTR(5,99) 84 stop 85 101 continue 86 do I=1,NA 87 do J=1,NRDF 88 IRDF(J,I)=0 89 end do 90 end do 91 do I=1,NAI 92 do J=1,NRDF 93 IRDFI(J,I)=0 94 end do 95 end do 96* Starting trajectory analysis 97 IEND=0 98 do while(IEND.eq.0) ! over configurations 99 50 call READCONF(IEND) 100 if(LCOMRDF1.or.LCOMRDF2)call GETCOM 101 if(IEND.ne.0)go to 200 102 if(ME.le.0)then 103 if(IPRINT.ge.6)write(*,'(a,f12.3,a,3f10.3)') 104 + ' Processing time',FULTIM 105 else 106 if(IPRINT.ge.6)write(*,'(a,f12.3,a,i5)') 107 + ' Processing time',FULTIM,' Ens ',MEE 108 if(ME.ne.MEE)then 109 IAN=IAN-1 110 go to 50 111 end if 112 end if 113 do I=1,NRDF ! over RDF groups 114 do K=1,IREP(I) ! over site pairs of the same group 115 IS=IRD(I,K) 116 JS=JRD(I,K) 117 ITYP=ITS(IS) 118 JTYP=ITS(JS) 119 NIS=NSPEC(ITYP) 120 NJS=NSPEC(JTYP) 121 do IM=1,NIS ! over first molecules 122 JBEG=1 123 if(IS.eq.JS)JBEG=IM+1 124 ISHIFT=ISADDR(ITYP)+IS-ISADR(ITYP)-NSITS(ITYP) 125 do JM=JBEG,NJS ! over second molecules 126 JSHIFT=ISADDR(JTYP)+JS-ISADR(JTYP)-NSITS(JTYP) 127 IST=ISHIFT+IM*NSITS(ITYP) 128 JST=JSHIFT+JM*NSITS(JTYP) 129 IMOL=NNUM(IST) 130 JMOL=NNUM(JST) 131 if(LCOMRDF1.and.LCOMRDF2)then 132 DX=X(IMOL)-X(JMOL) 133 DY=Y(IMOL)-Y(JMOL) 134 DZ=Z(IMOL)-Z(JMOL) 135 else if(LCOMRDF1.and..not.LCOMRDF2)then 136 DX=X(IMOL)-SX(JST) 137 DY=Y(IMOL)-SY(JST) 138 DZ=Z(IMOL)-SZ(JST) 139 else if(.not.LCOMRDF1.and.LCOMRDF2)then 140 DX=SX(IST)-X(JMOL) 141 DY=SY(IST)-Y(JMOL) 142 DZ=SZ(IST)-Z(JMOL) 143 else 144 DX=SX(IST)-SX(JST) 145 DY=SY(IST)-SY(JST) 146 DZ=SZ(IST)-SZ(JST) 147 end if 148 call PBC(DX,DY,DZ) 149 RR=sqrt(DX**2+DY**2+DZ**2) 150** write(*,*)I,K,IS,JS,IST,JST,IMOL,JMOL,RR 151* This distinguishes sites on the same molecule 152 if(IMOL.ne.JMOL)then 153 NR=NA*RR/RDFCUT+1 154 if(NR.gt.0.and.NR.le.NA)IRDF(I,NR)=IRDF(I,NR)+1 155 else 156 NR=NAI*(RR-RMI)/DRI+1 157** write(*,'(7i6,f8.4)')I,K,IS,JS,IST,JST,NR,RR 158 if(NR.gt.0.and.NR.le.NAI)IRDFI(I,NR)=IRDFI(I,NR)+1 159 end if 160 end do 161 end do 162 end do 163 end do 164 BOXS=BOXS+BOXL 165 BOYS=BOYS+BOYL 166 BOZS=BOZS+BOZL 167 end do 168 200 write(*,*)IAN,' configurations analysed' 169* RDF output 170 open(unit=16,file=FOUTRDF,status='unknown') 171 write(16,'(a)')' Radial distribution functions ' 172 write(16,*)IAN,' configurations analysed' 173* Averagre box sizes and volume 174 BOXS=BOXS/IAN 175 BOYS=BOYS/IAN 176 BOZS=BOZS/IAN 177 VOLS = BOXS*BOYS*BOZS 178 write(16,'(a,3f12.4)')' Average box sizes: ',BOXS,BOYS,BOZS 179 IFST=NA 180 do K=1,NRDF 181 CONU=0. 182 NPAIR=0 183 NPRI=0 ! for intramolecular. pairs 184 do KK=1,IREP(K) 185 IS=IRD(K,KK) 186 JS=JRD(K,KK) 187 ITYP=ITS(IS) 188 JTYP=ITS(JS) 189 NN1 = NSPEC(ITYP) 190 NN2 = NSPEC(JTYP) 191 if(ITYP.eq.JTYP)then 192 if(IS.eq.JS)then 193 NPAIR=NPAIR+NN1*(NN1-1)/2 194 else 195 NPAIR=NPAIR+NN1*(NN1-1) 196 end if 197 else 198 NPAIR=NPAIR+NN1*NN2 199 end if 200 if(ITYP.eq.JTYP.and.IS.ne.JS)NPRI=NPRI+1 201 end do 202 write(16,*) 203 FACT=VOLS/(NPAIR*1.d0*IAN) 204 if(LCOMRDF1.and.LCOMRDF2)then 205 write(16,'(5a,I12,a,I5)')' This mol. pair: ',NAME(ITS(IS)), 206 +' - ',NAME(ITS(JS)),' Npairs=',NPAIR,' Nint=',NPRI 207 else if(LCOMRDF1.and..not.LCOMRDF2)then 208 write(16,'(5a,I12,a,I5)')' This pair: ',NAME(ITS(IS)), 209 +' - ',NM(JS),' Npairs=',NPAIR,' Nint=',NPRI 210 else if(.not.LCOMRDF1.and.LCOMRDF2)then 211 write(16,'(5a,I12,a,I5)')' This pair: ',NM(IS), 212 +' - ',NAME(ITS(JS)),' Npairs=',NPAIR,' Nint=',NPRI 213 else 214 write(16,'(5a,I12,a,I5)')' This pair: ',NM(IS),' - ',NM(JS), 215 + ' Npairs=',NPAIR,' Nint=',NPRI 216 end if 217 write(16,'(a)') 218 + ' R RDF RCN1 RCN2 KBI(nm^3)' 219 write(16,*)'-------------------------------------------------' 220 CONU=0. 221 RKBI=0. 222 do I=1,NA 223 RDFV=IRDF(K,I) 224 RR=(I-0.5)*RDFCUT/NA 225 DRR=RDFCUT/NA 226 SH12=4*PI*DRR*(RR**2+DRR**2/12.d0) 227 CONU=CONU+RDFV*FACT/VOLS 228 CONU1=CONU*NSPEC(ITS(IS)) 229 CONU2=CONU*NSPEC(ITS(JS)) 230 GofR = RDFV*FACT/SH12 231 RKBI=RKBI+(GofR-1.)*SH12*0.001 232 if(CONU.gt.1.d-18)then 233 if(LCOMRDF1.and.LCOMRDF2)then 234 write(16,'(f8.3,3f12.5,2x,f12.5,3x,a4,2x,a4,a1,a4,I12)') 235 +RR,GofR,CONU1,CONU2,RKBI,'rdf:',NAME(ITS(IS)),'-', 236 +NAME(ITS(JS)),IRDF(K,I) 237 else if(LCOMRDF1.and..not.LCOMRDF2)then 238 write(16,'(f8.3,3f12.5,2x,f12.5,3x,a4,2x,a4,a1,a4,I12)') 239 +RR,GofR,CONU1,CONU2,RKBI,'rdf:',NAME(ITS(IS)),'-', 240 +NM(JS),IRDF(K,I) 241 else if(.not.LCOMRDF1.and.LCOMRDF2)then 242 write(16,'(f8.3,3f12.5,2x,f12.5,3x,a4,2x,a4,a1,a4,I12)') 243 +RR,GofR,CONU1,CONU2,RKBI,'rdf:',NM(IS),'-', 244 +NAME(ITS(JS)),IRDF(K,I) 245 else 246 write(16,'(f8.3,3f12.5,2x,f12.5,3x,a4,2x,a4,a1,a4,I12)') 247 +RR,GofR,CONU1,CONU2,RKBI,'rdf:',NM(IS),'-',NM(JS),IRDF(K,I) 248 end if 249 if(I.lt.IFST)IFST=I 250 end if 251 end do 252 if(NPRI.ge.1.and..not.(LCOMRDF1.and.LCOMRDF2))then 253 write(16,*) 254 write(16,*)' Intramolecular RDF ' 255 write(16,*)' ----------------------------------------- ' 256 write(16,'(a)') 257 + ' R RDF dist.distribution' 258 write(16,*)'-------------------------------------------------' 259 FACTI=1.d0/(NSPEC(ITS(IS))*NPRI*IAN) 260 do I=1,NAI 261 RR=RMI+(I-0.5)*DRI/NAI 262 DRR=DRI/NAI 263 SH12=4*PI*DRR*(RR**2+DRR**2/12.d0) 264 DIST=IRDFI(K,I)*FACTI/DRR 265 GofR=IRDFI(K,I)*FACT/SH12 266 if(DIST.gt.1.d-8)then 267 write(16,'(f8.3,2f13.5,11x,a8,2x,a4,a1,a4)') 268 + RR,GofR,DIST,'rdf_int:',NM(IS),'-',NM(JS) 269 end if 270 end do 271 end if 272 end do 273 write(16,*)' SUMMARY:' 274 write(16,'(a1,7x,50(1x,a4,a1,a4))') 275 + '#',(NM(IRD(K,1)),'-',NM(JRD(K,1)),K=1,NRDF) 276 do I=IFST,NA 277 RR=(I-0.5)*RDFCUT/NA 278 DRR=RDFCUT/NA 279 SH12=4*PI*DRR*(RR**2+DRR**2/12.d0) 280 do K=1,NRDF 281 NPAIR=0 282 do KK=1,IREP(K) 283 IS=IRD(K,KK) 284 JS=JRD(K,KK) 285 if(IS.eq.JS)then 286 NPAIR=NPAIR+NSPEC(ITS(IS))*(NSPEC(ITS(JS))-1)/2 287 else 288 NPAIR=NPAIR+NSPEC(ITS(IS))*NSPEC(ITS(JS)) 289 end if 290 end do 291 FACT=VOLS/(NPAIR*1.d0*IAN) 292 TMP(K)=IRDF(K,I)*FACT/SH12 293 end do 294 write(16,'(f8.3,50f10.5)')RR,(TMP(K),K=1,NRDF) 295 end do 296 close(16) 297 stop 298 end 299