1C 2C lgglib.f: functions used by Guoguang Lu's programs 3C Copyright (C) 1999 Guoguang Lu 4C 5C This library is free software: you can redistribute it and/or 6C modify it under the terms of the GNU Lesser General Public License 7C version 3, modified in accordance with the provisions of the 8C license to address the requirements of UK law. 9C 10C You should have received a copy of the modified GNU Lesser General 11C Public License along with this library. If not, copies may be 12C downloaded from http://www.ccp4.ac.uk/ccp4license.php 13C 14C This program is distributed in the hope that it will be useful, 15C but WITHOUT ANY WARRANTY; without even the implied warranty of 16C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 17C GNU Lesser General Public License for more details. 18C 19 function angle(v1,v2) 20c calculate the angle between two vectors v1, and v2 21c v1, v2 are the input 3*vectors 22c angle is the output angle between them (in degrees) 23 real*4 v1(3),v2(3) 24c 25 external acosd 26c write(7,*) 'v1,v2',v1,v2 27c write(7,*) 'vem v1,v2',vem(3,v1),vem(3,v2) 28c write(7,*) 'poimut',poimult(3,3,v1,v2) 29 arc = poimult(3,3,v1,v2)/(vem(3,v1)*vem(3,v2)) 30c write(7,*) 'arc=',arc 31 if (abs(arc).gt.1.) then 32 if (abs(arc)-1.gt.1e-5) write(7,*)'Warning arccosd > 1' 33 if (arc.gt.1.) arc = 1. 34 if (arc.lt.-1.) arc = -1. 35 end if 36 angle = acosd(arc) 37 end 38c 39c 40 SUBROUTINE ANGLZ(CELL,V,A) 41 DIMENSION CELL(6),V(3),A(3) 42 DO 10 I=1,3 4310 A(I)=CELL(I)*V(I) 44 END 45c 46c 47 SUBROUTINE ANTIARR(IM,IN,A1,A) 48c Transpose an array 49 DIMENSION A1(IM,IN),A(IN,IM) 50 DO 10 II=1,IM 51 DO 10 IJ=1,IN 5210 A(IJ,II)=A1(II,IJ) 53 END 54c 55c 56 SUBROUTINE ARRAD(IM,IN,A1,A2,OUT) 57c Add two arrays 58 DIMENSION A1(IM,IN),A2(IM,IN),OUT(IM,IN) 59 DO 10 II=1,IM 60 DO 10 IJ=1,IN 6110 OUT(II,IJ)=A1(II,IJ)+A2(II,IJ) 62 END 63c 64c 65 SUBROUTINE ARRGIVE(N,XYZ0,XYZ) 66C Copy an N-dimensional real array XYZ0 to another one XYZ. 67 REAL*4 XYZ(N),XYZ0(N) 68 DO I = 1, N 69 XYZ(I) = XYZ0(I) 70 END DO 71 END 72c 73c 74 SUBROUTINE ARRI(IM,IN,A,V,I) 75c Extract the i'th row of an array 76 INTEGER IM,IN 77 DIMENSION A(IM,IN) 78 DIMENSION V(IN) 79 DO 20 I1=1,IN 8020 V(I1)=A(I,I1) 81 RETURN 82 END 83c 84c 85 SUBROUTINE ARRJ(IM,IN,A,VJ,IJ) 86c Extract the j'th column of an array 87 DIMENSION A(IM,IN) 88 DIMENSION VJ(IM) 89 DO 20 I1=1,IM 9020 VJ(I1)=A(I1,IJ) 91 RETURN 92 END 93c 94c 95C------MULTIPLY A REAL ARRAY BY A CONSTANT 96 SUBROUTINE ARRMC(IM,IN,A1,C,A) 97 DIMENSION A1(IM,IN),A(IM,IN) 98 DO 10 I1=1,IM 99 DO 10 I2=1,IN 10010 A(I1,I2)=A1(I1,I2)*C 101 END 102c 103c 104C-------A STANDARD SUBPROGRAM TO MULTIPLY TWO ARRAYS 105C A1(IM1,IN1)*A2(IM2,IN2)=RSLT(IM1,IN2) 106 SUBROUTINE ARRMULT(IM1,IN1,IM2,IN2,A1,A2,RSLT,V1,V2) 107 INTEGER IM1,IM2,IN1,IN2 108 DIMENSION A1(IM1,IN1),A2(IM2,IN2),RSLT(IM1,IN2) 109 DIMENSION V1(IN1),V2(IM2) 110 IF (IN1.NE.IM2) WRITE(6,*) 'The two arrays cannot be multiplied' 111 DO 50 I=1,IM1 112 CALL ARRI(IM1,IN1,A1,V1,I) 113 DO 40 IJ=1,IN2 114 CALL ARRJ(IM2,IN2,A2,V2,IJ) 11540 RSLT(I,IJ)=POIMULT(IN1,IM2,V1,V2) 11650 CONTINUE 117 RETURN 118 END 119 120 SUBROUTINE ARRPS(IM,IN,A1,A2,OUT) 121 DIMENSION OUT(IM,IN),A1(IM,IN),A2(IM,IN) 122 DO 10 II=1,IM 123 DO 10 IJ=1,IN 12410 OUT(II,IJ)=A1(II,IJ)-A2(II,IJ) 125 END 126c 127c 128 subroutine arrrtoi(im,in,rmat,imat,err) 129c convert a matrix from real to integer 130c IM and IN are the dimensions. 131c rmat is input real IM*IN-dimension matrix 132c imat is output integer IM*IN-dimension matrix 133c err: if difference between real and integer values is bigger than 134c this parameter program will warn you 135 real*4 rmat(im,in) 136 integer imat(im,in) 137c 138 do j = 1, in 139 do i = 1, im 140 if (rmat(i,j).ge.0) then 141 imat(i,j) = int(rmat(i,j)+0.5) 142 else 143 imat(i,j) = int(rmat(i,j)-0.5) 144 end if 145 if (abs(rmat(i,j)-float(imat(i,j))).gt.err) then 146c WRITE(6,*) 'Warning: r-i > ',err 147c WRITE(6,*) 'i,j,real,integer',i,j,rmat(i,j),imat(i,j) 148 end if 149 end do 150 end do 151 end 152c 153c 154 subroutine arrvalue(n,array,value) 155c initialise an n-dimensional array 156 real array(n) 157 do i = 1, n 158 array(i) = value 159 end do 160 end 161c 162c 163 SUBROUTINE AVERG(M,N,XIN,AVE) 164C A routine to calculate the row averages of an M*N array XIN. 165C The mean M-dimensional array is AVE which average the N data in XIN. 166 DIMENSION XIN(M,N),AVE(M) 167 XN = N 168 DO I = 1, M 169 AVE(I) = 0. 170 DO J = 1, N 171 AVE(I) = AVE(I) + XIN(I,J) 172 END DO 173 AVE(I) = AVE(I)/ XN 174 END DO 175 END 176c 177c 178 function bondangle(xyz1,xyz2,xyz3) 179c subroutine to calculate a bond angle by giving 3 atoms. the output angle 180c will be vector 2-1 and 2-3 181c 1 3 182c \ / 183c \2/ 184c bondangle 185c xyz1,xyz2,xyz3 are the 3 atom positions and bondangle is shown in the fig 186c real xyz1(3),xyz2(3),xyz3(3) 187 real xyz1(3),xyz2(3),xyz3(3) 188 real b1(3),b2(3) 189 call arrps(3,1,xyz1,xyz2,b1) 190 call arrps(3,1,xyz3,xyz2,b2) 191 bondangle = angle(b1,b2) 192 end 193c 194c 195 function bonddihed(xyz1,xyz2,xyz3,xyz4) 196c subroutine to calculate a dihedral angle by giving 4 atoms. 197c the output angle will be between plane of 1-2-3 and 2-3-4 198c 1 4 1 199c \ / \ 200c \2---3/ \2---3\ bonddihed=180 in this case 201c bonddihed= 0 in this case \ 202c 4 203c xyz1,xyz2,xyz3 xyz4 are the 4 atom positions and bonddihed is 204c shown in the fig 205c real xyz1(3),xyz2(3),xyz3(3),xyz4(3) 206 real xyz1(3),xyz2(3),xyz3(3),xyz4(3) 207c real b1(3),b2(3) 208 real b321(3),b432(3) 209 real b23(3),b32(3),b34(3),b21(3) 210 real by(3),bdih(3),bd(2) 211 real view(3,3),views(3,3) 212 external acosd, asind 213c 214 call arrps(3,1,xyz3,xyz2,b23) 215 call arrps(3,1,xyz1,xyz2,b21) 216 call veccrsmlt(b23,b21,b321) 217 call arrps(3,1,xyz2,xyz3,b32) 218 call arrps(3,1,xyz4,xyz3,b34) 219 call veccrsmlt(b34,b32,b432) 220c call arrgive(3,b321,view(1,1)) 221c call arrgive(3,b23,view(1,3)) 222 call veccrsmlt(b32,b321,by) 223 call arrmc(3,1,b321,1./vem(3,b321),view(1,1)) 224 call arrmc(3,1,b32,1./vem(3,b32),view(1,3)) 225 call arrmc(3,1,by,1./vem(3,by),view(1,2)) 226 call antiarr(3,3,view,views) 227 call matmult(3,3,3,1,views,b432,bdih) 228 if (bdih(3).gt.1.0e-5) then 229 WRITE(6,*) 'Warning in dihedral',bdih 230 end if 231c 232 call arrmc(2,1,bdih,1./vem(2,bdih),bd) 233c 234c dump here 235c WRITE(6,*) 'views',views 236c WRITE(6,*) 'bdih',bdih 237c WRITE(6,*) 'bd',bd 238c 239 if (bd(1).ge.0. .and. bd(2).ge. 0.) then 240 bonddihed = acosd(bd(1)) 241 else if (bd(1).ge.0. .and. bd(2).lt.0) then 242 bonddihed = asind(bd(2)) 243 else if (bd(1).lt.0. .and. bd(2).ge.0) then 244 bonddihed = acosd(bd(1)) 245 else if (bd(1).lt.0. .and. bd(2).lt.0) then 246 bonddihed = -acosd(bd(1)) 247c else 248c stop ' what could be?' 249 end if 250 end 251c 252c 253 subroutine caseres(res1,res2) 254c A subroutine to change the case of amino acid e.g. PHE to Phe 255c res1 input residue name 256c res2 output residue name 257c 258 character*3 res1,res2,nam1(20),nam2(20) 259c... data statements. Separate declaration and init required for f2c 260 data nam1 / 261 1 'ALA','ARG','ASN','ASP','CYS','GLN','GLU','GLY','HIS','ILE', 262 2 'LEU','LYS','MET','PHE','PRO','SER','THR','TRP','TYR','VAL'/ 263 data nam2 / 264 1 'Ala','Arg','Asn','Asp','Cys','Gln','Glu','Gly','His','Ile', 265 2 'Leu','Lys','Met','Phe','Pro','Ser','Thr','Trp','Tyr','Val'/ 266c 267 do i = 1,20 268 if (res1.eq.nam1(i)) then 269 res2 = nam2(i) 270 return 271 end if 272 end do 273 res2 = res1 274 end 275c 276c 277 SUBROUTINE CLLZ(NSYM,SYM,CELL) 278 DIMENSION SYM(3,4,NSYM),CELL(6) 279 DO 20 IS=1,NSYM 280 DO 20 IC=1,3 28120 SYM(IC,4,IS)=SYM(IC,4,IS)*CELL(IC) 282 END 283c 284c 285 SUBROUTINE COMPARE ( NATM, XYZ1, XYZ2, A, T ) 286C DIMENSION XYZ1(3,NATM), XYZ2(3,NATM), A(3,3), T(3) 287C A subroutine to compare two structure without superposition. 288c MOL2 i.e. XYZ2(3*NATM) = A(3*3) * XYZ1(3*NATM) + T(3) 289C Where NATM represents the number of the atoms in each molecule 290c which will be superimposed. 291c XYZ1(3,NATM) represents the coordinates of MOL1 292c XYZ2(3,NATM) represents the coordinates of MOL2 293C A(3*3) represents the output matrix. 294c T(3*3) represents the output vector. 295c 296C NATM can not be larger than NATM0 (currently 50000) or the parameter NATM0 297C should be modified in this routine. 298c 299C COMMON/SHIFS/ SCALR,SCALT 300c DATA SCALR/1./,SCALT/1./ 301C Note: There should be a scale of the shifts. i.e. Shift= scale * shifts 302c Here SCALR is the rotation shiftscale 303c Here SCALT is the translation shiftscale 304c The initial and suggested SCALR is 1. 305c The initial and suggested SCALt is 1. 306C 307C COMMON/RMS/ RMS,SMEAN,NREF,NREF1 308C RMS is the final R.M.S between two molecules. 309C SMEAN is the mean distance. 310c NREF is the number of cycles of rotation refinement 311c NREF1 is not used by this routine. 312c They could be used to judge the SCALR and SCALS 313c 314C COMMON/IAT/ IAT1,IAT2,IAT3,IAT 315C DATA SCALR/1./,IAT/0/ 316c by Guoguang Lu 317C 01/31/1995 at Purdue 318C 319 COMMON/RMS/ RMS,SMEAN,NREF,NREF1 320c COMMON/SHIFS/ SCALR,SCALS 321c COMMON/IAT/ IAT1,IAT2,IAT3,IAT 322c DATA SCALR/1./,SCALS/1./,IAT/0/,IAT1/1/,IAT2/2/,IAT3/3/ 323 PARAMETER (NATM0 = 50000) 324c NATM0 is the largest number of atoms. 325 DIMENSION XYZ1(3,NATM), XYZ2(3,NATM), A(3,3), T(3) 326 DIMENSION X1(3,NATM0),X2(3,NATM0),DIS(3,NATM0) 327 DIMENSION CEN1(3),CEN2(3) 328 DIMENSION B1(3),B2(3),B3(3),T0(3) 329c 330c DEG=180./3.14159265359 331c 332C Number of atoms cannot be more than NATM0 or less than 3. 333c 334C 335 CALL RTMOV(NATM,XYZ1,A,T,DIS) 336 CALL ARRPS(3,NATM,DIS,XYZ2,DIS) 337C 338 ATM = NATM 339 SMEAN = 0 340 DO I=1, NATM 341 D = VEM(3,DIS(1,I)) 342 SMEAN = SMEAN + D 343 END DO 344 SMEAN = SMEAN /ATM 345C 346c SMEAN = SIGMA (DISTANCE) / N 347C i i 348 RMS = SQRT ( DOSQ(NATM*3,DIS) / ATM ) 349C 350c RMS = SQRT( SIGMA (DISTANCE^2) / N ) 351C i i 352c 353 WRITE(6,*) 354 WRITE(6,*) 'R.M.S.' 355 WRITE(6,*) ' natm' 356 WRITE(6,*) 'SQRT( SIGMA (Distance(i)^2)/natm ) = ',RMS 357 WRITE(6,*) ' i=1' 358C 359 WRITE(6,*) 360 WRITE(6,*) 'Mean Distance:' 361 WRITE(6,*) ' natm' 362 WRITE(6,*) 'SIGMA (Distance(i)) /natm = ',SMEAN 363 WRITE(6,*) ' i=1' 364C 365 WRITE(6,*) 366 WRITE(6,*) 'Mol1 is superposed to Mol2.' 367 WRITE(6,*) 'The matrix and the vector are:' 368 WRITE(6,*) 369 WRITE(6,1) (A(1,J),J=1,3),CEN1(1),T0(1) 370 WRITE(6,2) (A(2,J),J=1,3),CEN1(2),T0(2) 371 WRITE(6,1) (A(3,J),J=1,3),CEN1(3),T0(3) 3721 FORMAT(1X,' (',3F10.6,' ) ( ',f10.5,' ) (' 373 1 ,f10.5,' )') 3742 FORMAT(1X,' X2 = (',3F10.6,' ) * ( X1 -',f10.5,' ) + (' 375 1 ,f10.5,' )') 376C 377 CALL MATMULT(3,3,3,1,A,CEN1,B3) 378 CALL ARRPS(3,1,T0,B3,B3) 379 CALL ARRAD(3,1,CEN1,B3,T) 380C 381 WRITE(6,*) 382 WRITE(6,*) 383 WRITE(6,4) (A(1,J),J=1,3),T(1) 384 WRITE(6,5) (A(2,J),J=1,3),T(2) 385 WRITE(6,4) (A(3,J),J=1,3),T(3) 3864 FORMAT(1X,' (',3F10.6,' ) ( ) (',f10.5,' )') 3875 FORMAT(1X,' X2 = (',3F10.6,' ) * ( X1 ) + (',f10.5,' )') 388C 389 END 390c 391c 392 SUBROUTINE CRSTLARR(CELL,CRST,CRIS) 393 DIMENSION CRST(3,3),CRIS(3,3),BUF1(3),BUF2(3), 394 1IUF(3),CELL(6) 395 INTEGER IUF 396 external cosd, sind 397 AP=CELL(4) 398 BA=CELL(5) 399 GA=CELL(6) 400 CALL ARRMC(3,1,CRIS,0.,CRIS) 401 COASTAR=(COSD(BA)*COSD(GA)-COSD(AP))/(SIND(BA) 402 1*SIND(GA)) 403 SIASTAR=SQRT(1.-COASTAR*COASTAR) 404 CRIS(1,1)=1. 405 CRIS(1,2)=COSD(GA) 406 CRIS(1,3)=COSD(BA) 407 CRIS(2,2)=SIND(GA) 408 CRIS(2,3)=-SIND(BA)*COASTAR 409 CRIS(3,3)=SIND(BA)*SIASTAR 410 CALL ARRMC(3,3,CRIS,1.,CRST) 411 CALL IVSN(3,CRST,BUF1,BUF2,IUF,VAL,1E-6) 412 RETURN 413 END 414c 415c 416 SUBROUTINE CRTMOVE(NATM,XIN,A,T0,T,XOUT) 417C A routine to compute: XOUT = A * (XIN-T0) + T 418C where XIN is the input 3*NATM array 419C XOUT is the output 3*NATM array 420c A(3,3) is a 3*3 rotation matrix 421c T0(3) is a 3-dimensional translational vector. 422c T(3) is a 3-dimensional translational vector. 423C Note: XIN has been changed into XIN-T0 after this subroutine. 424c 89/11/06, BMC,UU,SE 425c DIMENSION XIN(3,NATM),XOUT(3,NATM),A(3,3),T(3) 426C 427 DIMENSION XIN(3,NATM),XOUT(3,NATM),A(3,3),T(3),T0(3) 428 DO I = 1, NATM 429 CALL ARRPS(3,1,XIN(1,I),T0,XIN(1,I)) 430 END DO 431 CALL MATMULT(3,3,3,NATM,A,XIN,XOUT) 432 DO I = 1, NATM 433 CALL ARRAD(3,1,XOUT(1,I),T,XOUT(1,I)) 434 END DO 435 END 436c 437c 438 SUBROUTINE LGG_CRYSTAL(CELL,CELLS,DEOR,ORTH,DEORS,ORTHS) 439C PJX - renamed from CRYSTAL to avoid clashes with similarly named 440C common block in mapmask, ncsmask and dm 441C 442C this is a routine which inputs a cell parameter, CELL and calculates 443c CELLS*6 --- cell dimensions in reciprocal space 444c DEOR3*3 --- Deorthogonalization matrix in real space 445c ORTH3*3 --- Orthogonalization matrix in recepical space. 446c DEORS3*3 --- Deorthogonalization matrix in reciprocal space 447c ORTHS3*3 --- Orthogonalization matrix in reciprocal space. 448c All the matrices have cell dimension information in them. 449c Guoguang 931118 450c 451 REAL DEOR(3,3),ORTH(3,3),CELL(6) 452 REAL DEORS(3,3),ORTHS(3,3),CELLS(6) 453 REAL BUFF(3,3) 454 DIMENSION B1(3),B2(3),IUF(3) 455 external cosd, sind 456 AP=CELL(4) 457 BA=CELL(5) 458 GA=CELL(6) 459 COASTAR=(COSD(BA)*COSD(GA)-COSD(AP))/(SIND(BA) 460 1*SIND(GA)) 461 SIASTAR=SQRT(1.-COASTAR*COASTAR) 462 CALL ARRVALUE(9,ORTH,0.) 463 ORTH(1,1)=CELL(1) 464 ORTH(1,2)=CELL(2)*COSD(GA) 465 ORTH(2,2)=CELL(2)*SIND(GA) 466 ORTH(1,3)=CELL(3)*COSD(BA) 467 ORTH(2,3)=-CELL(3)*SIND(BA)*COASTAR 468 ORTH(3,3)=CELL(3)*SIND(BA)*SIASTAR 469C 470 CALL ARRGIVE(9,ORTH,DEOR) 471 CALL IVSN(3,DEOR,B1,B2,IUF,DE,1.0E-6) 472 CALL ANTIARR(3,3,ORTH,DEORS) 473 CALL ANTIARR(3,3,DEOR,ORTHS) 474C 475 DO I = 1, 3 476 CELLS(I) = VEM(3,orths(1,I)) 477 END DO 478C 479 CELLS(4) = ANGLE(ORTHS(1,2),ORTHS(1,3)) 480 CELLS(5) = ANGLE(ORTHS(1,3),ORTHS(1,1)) 481 CELLS(6) = ANGLE(ORTHS(1,1),ORTHS(1,2)) 482 RETURN 483 END 484c 485c 486 SUBROUTINE CRYSTREC(CELL,CELLS,DEOR,ORTH,DEORS,ORTHS) 487C this is a routine which inputs a cell parameter, CELL and calculates 488c matrices for deorthogonalization and orthgonalization 489c The convention in this routine is x=a* y in plane of a*-b*, z is C direction 490c CELLS*6 --- cell dimensions in reciprocal space 491c DEOR3*3 --- Deorthogonalization matrix in real space 492c ORTH3*3 --- Orthogonalization matrix in reciprocal space. 493c DEORS3*3 --- Deorthogonalization matrix in reciprocal space 494c ORTHS3*3 --- Orthogonalization matrix in reciprocal space. 495c All the matrices have cell dimension information in them. 496c Guoguang 950914 497c 498 REAL DEOR(3,3),ORTH(3,3),CELL(6) 499 REAL DEORS(3,3),ORTHS(3,3),CELLS(6) 500 REAL BUFF(3,3),CELL1(6) 501 DIMENSION B1(3),B2(3),IUF(3) 502c Get a PDB convention 503 CALL LGG_CRYSTAL(CELL,CELLS,DEOR,ORTH,DEORS,ORTHS) 504 CALL LGG_CRYSTAL(CELLS,CELL1,DEORS,ORTHS,DEOR,ORTH) 505c 506 end 507c 508c 509c This is a subroutine to split one line into multiple lines by detecting 510c a separation character for example ; | and so on 511c the multiple lines will be written to unit iout 512c iout -- unit to be written to 513c il -- length of the line 514c line -- the line of characters to be split 515c sep -- "separation character" 516c iline -- output number of lines after splitting. 517c guoguang 940616 518 subroutine cutline(iout,line,sep,iline) 519 character*(*) line 520 character*1 sep 521 iline = 0 522 is = 1 523c WRITE(6,*) line 524c WRITE(6,*) 'sep ',sep 525 il = lenstr(line) 526 do it=il,1,-1 527 if (ichar(line(it:it)).le.31) then 528 line(it:it) = ' ' 529 else 530 goto 9 531 end if 532 end do 5339 il = it 534c WRITE(6,*) 'il',il,line(1:il) 535 rewind (iout) 536 if (il.eq.0) return 53710 ip = index(line(is:il),sep) 538 if (ip.eq.0) then 539 il1 = lenstr(line(is:il)) 540 if (il1.gt.0) then 541 do isps = is, is+il1-1 542 if (line(isps:isps).ne.' ') then 543 write(iout,'(a)') line(isps:is+il1-1) 544 iline = iline + 1 545 rewind (iout) 546 return 547 end if 548 end do 549 stop 'Strange' 550 end if 551 else 552 if (ip.gt.1) then 553 il1 = lenstr(line(is:ip-2+is)) 554 if (il1.gt.0) then 555 do isps = is, is+il1-1 556 if (line(isps:isps).ne.' ') then 557 write(iout,'(a)') line(isps:is+il1-1) 558 iline = iline + 1 559 is = ip + is 560 if (is.le.il) goto 10 561 rewind (10) 562 return 563 end if 564 end do 565 stop 'Strange' 566 end if 567 end if 568 is = ip + is 569 if (is.le.il) goto 10 570 rewind (10) 571 return 572 end if 573C stop 'CUTLINE' 574 end 575c 576c 577 function dedx2(x1,x2,x3,y1,y2,y3,abc) 578c this is a program which calculates dy/dx(x2) by fitting 579c a quadratic through 3 points. 580c y = a*x^2 + b*x + c 581c dy/dx = 2*a*x + b 582c Here dedv = dy/dx(x2) = 2*a*x2 + b 583c When there are 3 points x1,y1, x2,y2, x3,y3, the coefficients a,b,c 584c can be obtained from a linear equation and dy/dx(x2) can thus be 585c obtained. 586c 587c 588 dimension abc(3),y(3),arr(3,3) 589 dimension b1(3),b2(3),m(3) 590c 591 y(1) = y1 592 y(2) = y2 593 y(3) = y3 594 arr(1,1) = x1 * x1 595 arr(2,1) = x2 * x2 596 arr(3,1) = x3 * x3 597 arr(1,2) = x1 598 arr(2,2) = x2 599 arr(3,2) = x3 600 arr(1,3) = 1. 601 arr(2,3) = 1. 602 arr(3,3) = 1. 603 call ivsn(3,arr,b1,b2,m,de,1e-5) 604 call matmult(3,3,3,1,arr,y,abc) 605 dedx2 = 2.*abc(1)*x2 + abc(2) 606 return 607 end 608c 609c 610c a subroutine to convert denzo missetting angle to a matrix 611c definition [rot]=[rotx]*[roty]*[rotz] 612c 613 subroutine denmis(phi,rot) 614 real phi(3),rot(3,3) 615 external cosd, sind 616 sinx = sind(phi(1)) 617 cosx = cosd(phi(1)) 618 siny = sind(phi(2)) 619 cosy = cosd(phi(2)) 620 sinz = sind(phi(3)) 621 cosz = cosd(phi(3)) 622 rot(1,1) = cosy*cosz 623 rot(2,1) = -cosx*sinz+sinx*siny*cosz 624 rot(3,1) = sinx*sinz+cosx*siny*cosz 625 rot(1,2) = cosy*sinz 626 rot(2,2) = cosx*cosz+sinx*siny*sinz 627 rot(3,2) = -sinx*cosz+cosx*siny*sinz 628 rot(1,3) = -siny 629 rot(2,3) = sinx*cosy 630 rot(3,3) = cosx*cosy 631 end 632c 633c 634 function dihedral(x1,x2,x3,x4) 635c calculate the dihedral angle of four input positions 636c x1, x2, x3, x4 are the four input positions 637c dihedral is the output dihedral angle in degrees 638 real*4 x1(3),x2(3),x3(3),x4(3) 639 real*4 b1(3),b2(3),b3(3),b4(3) 640c 641 call arrps(3,1,x1,x2,b1) 642 call arrps(3,1,x3,x2,b2) 643 call veccrsmlt(b1,b2,b3) 644 call arrps(3,1,x2,x3,b1) 645 call arrps(3,1,x4,x3,b2) 646 call veccrsmlt(b1,b2,b4) 647 dihedral = angle(b3,b4) 648 end 649c 650c 651 FUNCTION DIST(XYZ1,XYZ2) 652c Subroutine to calculate distance between two positions 653 REAL xyz1(3),xyz2(3) 654 DX = XYZ1(1)-XYZ2(1) 655 DY = XYZ1(2)-XYZ2(2) 656 DZ = XYZ1(3)-XYZ2(3) 657 DIST = SQRT(DX*DX + DY*DY + DZ*DZ) 658 END 659c 660c 661 FUNCTION DOSQ(N,V) 662c DIMENSION V(N) 663C DOSQ = SIGMA ( V(i)^2 ) 664C i 665 DIMENSION V(N) 666 DOSQ = 0. 667 DO 10 I=1,N 66810 DOSQ = DOSQ + V(I)*V(I) 669 RETURN 670 END 671c 672c 673c------------------------------------------------------------ 674 subroutine down(txt,len) 675c 676c convert character string to lower case 677c 678 character*(*) txt 679 character*80 save 680 character*26 ualpha,lalpha 681 data lalpha /'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/ 682 data ualpha /'abcdefghijklmnopqrstuvwxyz'/ 683 684 do 9000 i=1,len 685 if((txt(i:i).ge.'A').and.(txt(i:i).le.'Z')) then 686 match = index(lalpha,txt(i:i)) 687 save(i:i) = ualpha(match:match) 688 else 689 save(i:i) = txt(i:i) 690 end if 691c end do 6929000 continue 693 694 txt = save 695 return 696 end 697c 698c 699 SUBROUTINE DRVRBD(ITH,ROH,A) 700C This is a subroutine to calculate the matrix of deriv(E)/deriv(theta i) 701c where E is a 3*3 rotation matrix and theta i (i=1,2,3) is Eulerian Angle 702c in ROH system. 703C DIMENSION A(3,3),ROH(3) 704c parameter description: 705c ITH: when ITH = 1 706c the output matrix is deriv(E) /deriv(theta1) 707c when ITH = 2 708c the output matrix is deriv(E) /deriv(theta2) 709c when ITH = 3 710c the output matrix is deriv(E) /deriv(theta3) 711c ROH are the three Eulerian angles 712c A is the output deriv matrix. i.e. A = deriv(E) /deriv)(theta i) 713c written by Guoguang Lu 714c 715 DIMENSION A(3,3),ROH(3) 716 DIMENSION A1(3,3) 717 ARC=3.14159265359/180.0 718 SIPS = SIN(ROH(1)*ARC) 719 COPS = COS(ROH(1)*ARC) 720 SITH = SIN(ROH(2)*ARC) 721 COTH = COS(ROH(2)*ARC) 722 SIPH = SIN(ROH(3)*ARC) 723 COPH = COS(ROH(3)*ARC) 724c 725c Rossmann & Blow Matrix 726c 727c E(1,1) =-SIPS*COTH*SIPH + COPS*COPH 728c E(1,2) =-SIPS*COTH*COPH - COPS*SIPH 729c E(1,3) = SIPS*SITH 730c E(2,1) = COPS*COTH*SIPH + SIPS*COPH 731c E(2,2) = COPS*COTH*COPH - SIPS*SIPH 732c E(2,3) = - COPS*SITH 733c E(3,1) = SITH*SIPH 734c E(3,2) = SITH*COPH 735c E(3,3) = COTH 736 737 738 IF (ITH.EQ.1) THEN 739 A(1,1) =-coPS*COTH*SIPH - siPS*COPH 740 A(1,2) =-coPS*COTH*COPH + siPS*SIPH 741 A(1,3) = coPS*SITH 742 A(2,1) =-siPS*COTH*SIPH + coPS*COPH 743 A(2,2) =-siPS*COTH*COPH - coPS*SIPH 744 A(2,3) = + siPS*SITH 745 A(3,1) = 0. 746 A(3,2) = 0. 747 A(3,3) = 0. 748 ELSE IF (ITH.EQ.2) THEN 749 A(1,1) =+SIPS*siTH*SIPH 750 A(1,2) =+SIPS*siTH*COPH 751 A(1,3) = SIPS*coTH 752 A(2,1) =-COPS*siTH*SIPH 753 A(2,2) =-COPS*siTH*COPH 754 A(2,3) = - COPS*coTH 755 A(3,1) = coTH*SIPH 756 A(3,2) = coTH*COPH 757 A(3,3) = -siTH 758 ELSE IF (ITH.EQ.3) THEN 759 A(1,1) =-SIPS*COTH*coPH - COPS*siPH 760 A(1,2) =+SIPS*COTH*siPH - COPS*coPH 761 A(1,3) = 0. 762 A(2,1) = COPS*COTH*coPH - SIPS*siPH 763 A(2,2) =-COPS*COTH*siPH - SIPS*coPH 764 A(2,3) = 0. 765 A(3,1) = SITH*coPH 766 A(3,2) = -SITH*siPH 767 A(3,3) = 0. 768 ELSE 769 STOP 'invalid parameter ITH ' 770 END IF 771C SINCE ROH IS IN DEGREES, THE MATRIX HAS TO BE MULTIPLIED BY ARC 772 CALL ARRMC(3,3,A,ARC,A1) 773 call arrgive(9,a1,a) 774 END 775c 776c 777 SUBROUTINE DRVROHTH(ITH,ROH,A) 778C This is a subroutine to calculate the matrix of deriv(E)/deriv(theta i) 779c where E is a 3*3 rotation matrix and theta i (i=1,2,3) is Eulerian Angle 780c in ROH system. 781C DIMENSION A(3,3),ROH(3) 782c parameter description: 783c ITH: when ITH = 1 784c the output matrix is deriv(E) /deriv(theta1) 785c when ITH = 2 786c the output matrix is deriv(E) /deriv(theta2) 787c when ITH = 3 788c the output matrix is deriv(E) /deriv(theta3) 789c ROH are the three Eulerian angles 790c A is the output deriv matrix. i.e. A = deriv(E) /deriv)(theta i) 791c written by Guoguang Lu 792c 793 DIMENSION A(3,3),ROH(3) 794 ARC=3.14159265359/180.0 795 SIPS = SIN(ROH(1)*ARC) 796 COPS = COS(ROH(1)*ARC) 797 SITH = SIN(ROH(2)*ARC) 798 COTH = COS(ROH(2)*ARC) 799 SIPH = SIN(ROH(3)*ARC) 800 COPH = COS(ROH(3)*ARC) 801C 802C 803C * ROH - MATRIX 804C 805C 110 E(1,1) = COPS*COPH - SIPS*SIPH*SITH 806C E(1,2) =-SIPS*COTH 807C E(1,3) = COPS*SIPH + SIPS*SITH*COPH 808C E(2,1) = SIPS*COPH + COPS*SITH*SIPH 809C E(2,2) = COPS*COTH 810C E(2,3) = SIPS*SIPH - COPS*SITH*COPH 811C E(3,1) =-COTH*SIPH 812C E(3,2) = SITH 813C E(3,3) = COTH*COPH 814 815 IF (ITH.EQ.1) THEN 816 A(1,1) = -SIPS*COPH - COPS*SIPH*SITH 817 A(1,2) = -COPS*COTH 818 A(1,3) = -SIPS*SIPH + COPS*SITH*COPH 819 A(2,1) = COPS*COPH - SIPS*SITH*SIPH 820 A(2,2) = -SIPS*COTH 821 A(2,3) = COPS*SIPH + SIPS*SITH*COPH 822 A(3,1) = 0. 823 A(3,2) = 0. 824 A(3,3) = 0. 825 ELSE IF (ITH.EQ.2) THEN 826 A(1,1) = - SIPS*SIPH*COTH 827 A(1,2) = SIPS*SITH 828 A(1,3) = SIPS*COTH*COPH 829 A(2,1) = COPS*COTH*SIPH 830 A(2,2) = -COPS*SITH 831 A(2,3) = - COPS*COTH*COPH 832 A(3,1) = SITH*SIPH 833 A(3,2) = COTH 834 A(3,3) = -SITH*COPH 835 ELSE IF (ITH.EQ.3) THEN 836 A(1,1) = -COPS*SIPH - SIPS*COPH*SITH 837 A(1,2) = 0. 838 A(1,3) = COPS*COPH - SIPS*SITH*SIPH 839 A(2,1) =- SIPS*SIPH + COPS*SITH*COPH 840 A(2,2) = 0. 841 A(2,3) = SIPS*COPH + COPS*SITH*SIPH 842 A(3,1) =-COTH*COPH 843 A(3,2) = 0. 844 A(3,3) = - COTH*SIPH 845 ELSE 846 STOP 'invalid parameter ITH ' 847 END IF 848 END 849c 850c 851 SUBROUTINE DRVROHTHARC(ITH,ROH,A) 852C This is a subroutine to calculate the matrix of deriv(E)/deriv(theta i) 853c where E is a 3*3 rotation matrix and theta i (i=1,2,3) is Eulerian Angle 854c in ROH system. 855C DIMENSION A(3,3),ROH(3) 856c parameter description: 857c ITH: when ITH = 1 858c the output matrix is deriv(E) /deriv(theta1) 859c when ITH = 2 860c the output matrix is deriv(E) /deriv(theta2) 861c when ITH = 3 862c the output matrix is deriv(E) /deriv(theta3) 863c ROH are the three Eulerian angles 864c A is the output deriv matrix. i.e. A = deriv(E) /deriv)(theta i) 865c written by Guoguang Lu 866c 867 DIMENSION A(3,3),ROH(3) 868C ARC=3.14159265359/180.0 869C ARC= 1.0 870 SIPS = SIN(ROH(1)) 871 COPS = COS(ROH(1)) 872 SITH = SIN(ROH(2)) 873 COTH = COS(ROH(2)) 874 SIPH = SIN(ROH(3)) 875 COPH = COS(ROH(3)) 876C 877C 878C * ROH - MATRIX 879C 880C 110 E(1,1) = COPS*COPH - SIPS*SIPH*SITH 881C E(1,2) =-SIPS*COTH 882C E(1,3) = COPS*SIPH + SIPS*SITH*COPH 883C E(2,1) = SIPS*COPH + COPS*SITH*SIPH 884C E(2,2) = COPS*COTH 885C E(2,3) = SIPS*SIPH - COPS*SITH*COPH 886C E(3,1) =-COTH*SIPH 887C E(3,2) = SITH 888C E(3,3) = COTH*COPH 889 890 IF (ITH.EQ.1) THEN 891 A(1,1) = -SIPS*COPH - COPS*SIPH*SITH 892 A(1,2) = -COPS*COTH 893 A(1,3) = -SIPS*SIPH + COPS*SITH*COPH 894 A(2,1) = COPS*COPH - SIPS*SITH*SIPH 895 A(2,2) = -SIPS*COTH 896 A(2,3) = COPS*SIPH + SIPS*SITH*COPH 897 A(3,1) = 0. 898 A(3,2) = 0. 899 A(3,3) = 0. 900 ELSE IF (ITH.EQ.2) THEN 901 A(1,1) = - SIPS*SIPH*COTH 902 A(1,2) = SIPS*SITH 903 A(1,3) = SIPS*COTH*COPH 904 A(2,1) = COPS*COTH*SIPH 905 A(2,2) = -COPS*SITH 906 A(2,3) = - COPS*COTH*COPH 907 A(3,1) = SITH*SIPH 908 A(3,2) = COTH 909 A(3,3) = -SITH*COPH 910 ELSE IF (ITH.EQ.3) THEN 911 A(1,1) = -COPS*SIPH - SIPS*COPH*SITH 912 A(1,2) = 0. 913 A(1,3) = COPS*COPH - SIPS*SITH*SIPH 914 A(2,1) =- SIPS*SIPH + COPS*SITH*COPH 915 A(2,2) = 0. 916 A(2,3) = SIPS*COPH + COPS*SITH*SIPH 917 A(3,1) =-COTH*COPH 918 A(3,2) = 0. 919 A(3,3) = - COTH*SIPH 920 ELSE 921 STOP 'invalid parameter ITH ' 922 END IF 923 END 924c 925c 926 SUBROUTINE DRVROHTHD(ITH,ROH,A) 927C This is a subroutine to calculate the matrix of deriv(E)/deriv(theta i) 928c where E is a 3*3 rotation matrix and theta i (i=1,2,3) is Eulerian Angle 929c in ROH system. 930C DIMENSION A(3,3),ROH(3) 931c parameter description: 932c ITH: when ITH = 1 933c the output matrix is deriv(E) /deriv(theta1) 934c when ITH = 2 935c the output matrix is deriv(E) /deriv(theta2) 936c when ITH = 3 937c the output matrix is deriv(E) /deriv(theta3) 938c ROH are the three Eulerian angles 939c A is the output deriv matrix. i.e. A = deriv(E) /deriv)(theta i) 940c written by Guoguang Lu 941c 942 DIMENSION A(3,3),ROH(3) 943 DIMENSION A1(3,3) 944 ARC=3.14159265359/180.0 945 SIPS = SIN(ROH(1)*ARC) 946 COPS = COS(ROH(1)*ARC) 947 SITH = SIN(ROH(2)*ARC) 948 COTH = COS(ROH(2)*ARC) 949 SIPH = SIN(ROH(3)*ARC) 950 COPH = COS(ROH(3)*ARC) 951C 952C 953C * ROH - MATRIX 954C 955C 110 E(1,1) = COPS*COPH - SIPS*SIPH*SITH 956C E(1,2) =-SIPS*COTH 957C E(1,3) = COPS*SIPH + SIPS*SITH*COPH 958C E(2,1) = SIPS*COPH + COPS*SITH*SIPH 959C E(2,2) = COPS*COTH 960C E(2,3) = SIPS*SIPH - COPS*SITH*COPH 961C E(3,1) =-COTH*SIPH 962C E(3,2) = SITH 963C E(3,3) = COTH*COPH 964 965 IF (ITH.EQ.1) THEN 966 A(1,1) = -SIPS*COPH - COPS*SIPH*SITH 967 A(1,2) = -COPS*COTH 968 A(1,3) = -SIPS*SIPH + COPS*SITH*COPH 969 A(2,1) = COPS*COPH - SIPS*SITH*SIPH 970 A(2,2) = -SIPS*COTH 971 A(2,3) = COPS*SIPH + SIPS*SITH*COPH 972 A(3,1) = 0. 973 A(3,2) = 0. 974 A(3,3) = 0. 975 ELSE IF (ITH.EQ.2) THEN 976 A(1,1) = - SIPS*SIPH*COTH 977 A(1,2) = SIPS*SITH 978 A(1,3) = SIPS*COTH*COPH 979 A(2,1) = COPS*COTH*SIPH 980 A(2,2) = -COPS*SITH 981 A(2,3) = - COPS*COTH*COPH 982 A(3,1) = SITH*SIPH 983 A(3,2) = COTH 984 A(3,3) = -SITH*COPH 985 ELSE IF (ITH.EQ.3) THEN 986 A(1,1) = -COPS*SIPH - SIPS*COPH*SITH 987 A(1,2) = 0. 988 A(1,3) = COPS*COPH - SIPS*SITH*SIPH 989 A(2,1) =- SIPS*SIPH + COPS*SITH*COPH 990 A(2,2) = 0. 991 A(2,3) = SIPS*COPH + COPS*SITH*SIPH 992 A(3,1) =-COTH*COPH 993 A(3,2) = 0. 994 A(3,3) = - COTH*SIPH 995 ELSE 996 STOP 'invalid parameter ITH ' 997 END IF 998C SINCE ROH IS IN DEGREES, THE MATRIX HAS TO MULTIPLIED BY ARC 999 CALL ARRMC(3,3,A,ARC,A1) 1000 CALL arrgive(9,A1,A1) 1001 END 1002c 1003c 1004 function dstll2(x1,x2,y1,y2) 1005c a routine to calculate the line-line distance. 1006c input: 1007c x1 and x2 are two positions on line X 1008c y1 and y2 are two positions on line Y 1009c output 1010c dstll1 is the distance between x and y 1011 real x1(3),x2(3) 1012 real y1(3),y2(3) 1013 real b1(3),b2(3),b3(3),b4(3),b5(3),b6(3) 1014c 1015 call arrps(3,1,x2,x1,b1) 1016 call arrps(3,1,y2,y1,b2) 1017 call veccrsmlt(b1,b2,b3) 1018 if (vem(3,b3).ge.1e-6) then 1019 dstll2 = dstps1(b3,x1,y1) 1020 else 1021 call arrps(3,1,y1,x1,b5) 1022 call arrmc(3,1,b1,1./vem(3,b1),b4) 1023 d1 = poimult(3,3,b5,b4) 1024 d2 = vem(3,b5) 1025 dstll2 = sqrt(d2*d2-d1*d1) 1026 end if 1027 end 1028c 1029c 1030 FUNCTION DSTPL1(VL,XYZ0,XYZ,VPL) 1031C-----------DISTANCE FROM A POINT TO LINE 1032c FUNCTION DSTPL1(VL,XYZ0,XYZ,VPL) 1033c VL -3- is the input direction vector of the line 1034c XYZ0 -3- is the input position of a point on the line 1035c XYZ -3- is the input coordinates of the point. 1036c VPL -3- is the output vector from the point to the line which is 1037c perpendicular to the direction vector of the line. 1038c DISTPL1 is the output distance from the point to the line i.e. /VPL(3)/ 1039C so the equation of the line should be 1040C x-XYZ0(1) y-XYZ0(2) z-XYZ0(3) 1041C ----------- = ----------- = ----------- 1042C VL(1) VL(2) VL(3) 1043C where x,y,z is any position on the line 1044c 1045 1046 DIMENSION VL(3),XYZ0(3),XYZ(3),VPL(3) 1047 DIMENSION B(3),BUF(3) 1048 CALL ARRPS(3,1,XYZ,XYZ0,BUF) 1049 T=POIMULT(3,3,BUF,VL)/POIMULT(3,3,VL,VL) 1050 CALL ARRMC(3,1,VL,T,BUF) 1051 CALL ARRAD(3,1,BUF,XYZ0,B) 1052 CALL ARRPS(3,1,XYZ,B,VPL) 1053 DSTPL1=VEM(3,VPL) 1054 RETURN 1055 END 1056c 1057c 1058 FUNCTION DSTPL2(VL,XYZ0,XYZ,VP0,VPL) 1059C-----------DISTANCE FROM A POINT TO LINE 1060c FUNCTION DSTPL1(VL,XYZ0,XYZ,VP0,VPL) 1061c VL -3- is the input directional vector of the line 1062c XYZ0 -3- is the a input position which belong to the line 1063c XYZ -3- is the input coordinate of the point. 1064c VP0 -3- is the output coordinate whis is the image of the point in the 1065c line 1066C VPL -3- is the vector from VP0 to XYZ 1067c DISTPL1 is the output distance from the point to the line i.e. /VPL(3)/ 1068C so the equation of the line should be 1069C x-XYZ0(1) y-XYZ0(2) z-XYZ0(3) 1070C ----------- = ----------- = ----------- 1071C VL(1) VL(2) VL(3) 1072C where the x,y,z is any position of the line 1073c 1074 1075 DIMENSION VL(3),XYZ0(3),XYZ(3),VP0(3) 1076 DIMENSION VPL(3),BUF(3) 1077 CALL ARRPS(3,1,XYZ,XYZ0,BUF) 1078 DVL = VEM(3,VL) 1079 T= POIMULT(3,3,BUF,VL)/DVL*DVL 1080 CALL ARRMC(3,1,VL,T,BUF) 1081 CALL ARRAD(3,1,BUF,XYZ0,VP0) 1082 CALL ARRPS(3,1,XYZ,VP0,VPL) 1083 DSTPL2=VEM(3,VPL) 1084 RETURN 1085 END 1086c 1087c 1088 FUNCTION DSTPLROTR(A,T,XYZ) 1089C-----------DISTANCE FROM A POINT TO A ROTATION axis defined by ROT,T 1090c FUNCTION DSTPL1(VL,XYZ0,XYZ,VPL) 1091c A 3*3 is the rotation matrix which defines the axis 1092c T 3 is the translation vector which defines the axis 1093c When there is a non-crystallographic rotation and translation 1094c x2 = A*x1 + t 1095c where A(3,3) is the 3*3 matrix 1096c t(3) is translation vector. 1097c XYZ -3- is the input coordinate of the point. 1098c DISTPLROTR is the output distance from the point to the line i.e. /VPL(3)/ 1099C so the equation of the line should be 1100C where the x,y,z is any position of the line 1101c Guoguang 970921 1102c 1103 real a(3,3),t(3),xyz(3),POL(3),xyz0(3),vpl(3) 1104 real B(3),BUF(3),vl(3),vl0(3) 1105 real b1(3),b2(3),x1(3),x2(3) 1106 external sind 1107c 1108 call arrvalue(3,xyz0,0.) 1109 call mtovec(a,vl,vkapa) 1110 vlm = vem(3,vl) 1111 if (vlm.lt.1.e-4) then 1112 WRITE(6,*) 'a',a 1113 WRITE(6,*) 'vlm',vlm 1114 stop 'something strange in dstplrotr' 1115 end if 1116 if (abs(vlm-1.).gt.1.e-5) then 1117 call arrmc(3,1,vl,1./vlm,vl0) 1118 call arrgive(3,vl0,vl) 1119 end if 1120 if (abs(vkapa).lt.1e-4) then 1121 go = vem(3,t) 1122 dstplrotr=9999. 1123 else 1124 go = poimult(3,3,t,vl) 1125 end if 1126 call arrmc(3,1,vl,go,b1) 1127 call rtmov(3,xyz,a,t,x1) 1128 call arrps(3,1,x1,b1,x2) 1129 dstplrotr = dist(x2,xyz)/(2.*sind(vkapa/2.)) 1130 end 1131c 1132c 1133 FUNCTION DSTPS1(VS,XYZ0,XYZ) 1134C------DISTANCE FROM A POINT TO A SECTION 1135c FUNCTION DSTPS1(VL,XYZ0,XYZ,VPL) 1136c VL -3- is the input directional vector of the plane section 1137c XYZ0 -3- is the a input position which belong to the section 1138c XYZ -3- is the input coordinate of the point. 1139c DISTPL1 is the output distance from the point to the section 1140C so the equation of the section should be 1141C 1142c VS(1) * ( x - XYZ0(1) ) + VS(2) * ( y - XYZ0(2) ) + VS(3) ( z - XYZ0(3) ) = 0 1143C 1144C where the x,y,z is any position of the section 1145c 1146 DIMENSION VS(3),XYZ0(3),XYZ(3),BUF(3) 1147 CALL ARRPS(3,1,XYZ,XYZ0,BUF) 1148 DSTPS1=POIMULT(3,3,VS,BUF)/VEM(3,VS) 1149 RETURN 1150 END 1151c 1152c 1153c---------------------------------------------------------- 1154 subroutine elb(txt,len) 1155c 1156c eliminate leading blanks from a character string 1157c 1158 character*(*) txt 1159 character*80 save 1160 1161 do 9000 i=1,len 1162 if(txt(i:i).ne.' ') then 1163 nonb = i 1164 go to 100 1165 end if 11669000 continue 1167 return 1168100 continue 1169 save = txt(nonb:len) 1170 txt = save 1171 return 1172 end 1173c 1174c 1175 SUBROUTINE ELIZE(N,E) 1176 DIMENSION E(N,N) 1177 DO 10 I=1,N 1178 E(I,I)=1. 1179c WRITE(6,*) I,E 1180 DO 5 IJ=1,N 1181 IF (I.NE.IJ) E(I,IJ)=0. 11825 CONTINUE 118310 CONTINUE 1184 RETURN 1185 END 1186c 1187c 1188 function error2d(xy1,xy2) 1189 real xy1(2),xy2(2) 1190 e1 = xy2(1)-xy1(1) 1191 e2 = xy2(2)-xy1(2) 1192 error2d = sqrt(e1*e1+e2*e2) 1193 end 1194c 1195c 1196c read in a file and copy to another file 1197c nin --- unit of input file 1198c nout -- unit of output file 1199c file -- name of the input file 1200 subroutine filein(nin,nout,file) 1201 character*(*) file 1202 character*132 line 1203 call ccpdpn(nin,file,'readonly','F',0,0) 12045 read(nin,'(a)',end=10) line 1205 il = lenstr(line) 1206 write(nout,'(a)') line(1:il) 1207 goto 5 120810 continue 1209 end 1210c 1211c 1212 SUBROUTINE FMATIN(N,X) 1213 REAL X(N) 1214 DO I = 1, N 1215 CALL FRCINSIDE(X(i)) 1216 END DO 1217 END 1218c 1219c 1220 SUBROUTINE FRAC2ANG(XYZ) 1221C A subroutine to convert fractional coordinates xyz into orthogonal 1222c coordinates in Angstrom units. 1223C 1224C SUBROUTINE FRAC2ANG(XYZ) 1225C COMMON/CELL/ CELL,DEOR,ORTH 1226C DIMENSION CELL(6),DEOR(3,3),ORTH(3,3),XYZ(3) 1227C XYZ is an input fractional coordinate 3-dim array. 1228c CELL: 6-dim array to save the cell constants 1229c DEOR is 3*3 array, representing the deorthogonalizing array. 1230c ORTH is 3*3 array, representing the orthogonalizing array. 1231c 1232 COMMON/CELL/ CELL,DEOR,ORTH 1233 DIMENSION CELL(6),DEOR(3,3),ORTH(3,3),XYZ(3) 1234 DIMENSION B1(3),B2(3),B3(3) 1235 B3(1) = XYZ(1) * CELL(1) 1236 B3(2) = XYZ(2) * CELL(2) 1237 B3(3) = XYZ(3) * CELL(3) 1238C 1239 CALL matmult(3,3,3,1,ORTH,B3,XYZ,B1) 1240 END 1241c 1242c 1243 SUBROUTINE FRCINSIDE(X) 124410 CONTINUE 1245 IF (X.GE.1.) THEN 1246 X = MOD(X,1.) 1247 ELSE IF (X.LT.0.) THEN 1248 X = MOD(X,1.) + 1. 1249 END IF 1250 END 1251c 1252c 1253 SUBROUTINE FRCTOANG(XYZ) 1254C A subroutine to convert fractional coordinates xyz into orthogonal 1255c coordinates in Angstrom units. 1256C 1257C SUBROUTINE FRAC2ANG(XYZ) 1258C COMMON/CELL/ CELL,DEOR,ORTH 1259C DIMENSION CELL(6),DEOR(3,3),ORTH(3,3),XYZ(3) 1260C XYZ is an input fractional coordinate 3-dim array. 1261c CELL: 6-dim array to save the cell constants 1262c DEOR is 3*3 array, representing the deorthogonalizing array. 1263c ORTH is 3*3 array, representing the orthogonalizing array. 1264c obs: different from frac2ang, the deor and orth must be from 1265c subroutine crystal 1266c 1267 COMMON/CELL/ CELL,DEOR,ORTH 1268 DIMENSION CELL(6),DEOR(3,3),ORTH(3,3),XYZ(3) 1269 DIMENSION B1(3),B2(3),B3(3) 1270c B3(1) = XYZ(1) * CELL(1) 1271c B3(2) = XYZ(2) * CELL(2) 1272c B3(3) = XYZ(3) * CELL(3) 1273C 1274 B3(1) = XYZ(1) 1275 B3(2) = XYZ(2) 1276 B3(3) = XYZ(3) 1277 CALL matmult(3,3,3,1,ORTH,B3,XYZ,B1) 1278 END 1279 character*80 function getnam(filnam) 1280 character*(*) filnam 1281c character*80 cmd 1282c call spstrunct(filnam) 1283c il = lenstr(filnam) 1284c write(cmd,'(3a)') 'printenv ',filnam(1:il),' > GETNAM' 1285cc call system(cmd) 1286c call ccpdpn(31,'GETNAM','old','F',0,0) 1287c read(31,'(a)') getnam 1288c call spstrunct(getnam) 1289c close (31) 1290 getnam = ' ' 1291 call ugtenv(filnam,getnam) 1292 end 1293c 1294c 1295 SUBROUTINE GETPDB(X,ATOM,RES,SEQ,NATM,FILNAM) 1296C A subroutine to read the Protein data bank format coordinate file 1297c X 3*21000 arrays, output xyz coordinates 1298c ATOM 21000 character*4 arrays, output atom names 1299c RES 21000 character*4 arrays, output residue number e.g. A101 A103... 1300C If residue number is like A 1 then it will become A1 1301c SEQ 21000 character*4 arrays, output peptide name. e.g. PHE TRY .... 1302C NATM is the number of atoms in the coordinate file 1303c FILNAT is the file name of the coordinate. 1304c 1305 DIMENSION X(3,21000) 1306 CHARACTER*4 ATOM(21000),RES(21000),SEQ(21000),RES2 1307 CHARACTER FILNAM*80,HEAD*6,RES1*5 1308C DATA FILNAM/'DIAMOND'/ 1309C 1310 CALL SPSTRUNCT(FILNAM) 1311 call ccpdpn(3,FILNAM,'READONLY','F',0,0) 1312C 1313 NATM = 1 13145 READ(3,10,ERR=5,END=20) HEAD, ATOM(NATM), SEQ(NATM), RES1 1315 1 , (X(I,NATM),I=1,3) 131610 FORMAT(A6,7X,2A4,A5,4X,3F8.3) 1317 IF (HEAD.NE.'ATOM '.AND.HEAD.NE.'HEDATM') GOTO 5 1318C 1319 IF (RES1(1:1).NE.' '.AND.RES1(2:2).EQ.' ') THEN 1320 RES2(1:1) = RES1(1:1) 1321 RES2(2:4) = RES1(3:5) 1322 ELSE 1323 RES2 = RES1(2:5) 1324 END IF 1325C 1326 call spstrunct(atom(natm)) 1327 call spstrunct(seq(natm)) 1328 call spstrunct(RES2) 1329 res(natm)=res2 1330C 133112 CONTINUE 1332C 1333 NATM = NATM + 1 1334 IF (NATM.GT.21000) STOP 'Atoms can not be more than 21000.' 1335 GOTO 5 133620 CONTINUE 1337 NATM = NATM - 1 1338 CLOSE (3) 1339C 1340 END 1341c 1342c 1343 SUBROUTINE GETPDB1(X,ATOM,SEQ,RES,BFAC,OCC,NATM,FILNAM) 1344C A subroutine to read the Protein data bank format coordinate file 1345c X 3*15000 arrays, output xyz coordinates 1346c ATOM 15000 character*4 arrays, output atom names 1347c RES 15000 character*4 arrays, output residue number e.g. A101 A103... 1348C If residue number is like A 1 then it will become A1 1349c SEQ 15000 character*4 arrays, output peptide name. e.g. PHE TRY .... 1350C NATM is the number of atoms in the coordinate file 1351c FILNAT is the file name of the coordinate. 1352c 1353c in this subroutine, hydrogen atoms were excluded. 1354c 1355 PARAMETER (MAXATM = 25000) 1356 DIMENSION X(3,MAXATM) 1357 real bfac(MAXATM),occ(MAXATM) 1358 CHARACTER*4 ATOM(MAXATM),RES(MAXATM),SEQ(MAXATM),RES2 1359 CHARACTER FILNAM*80,HEAD*6,RES1*5 1360C DATA FILNAM/'DIAMOND'/ 1361C 1362 CALL SPSTRUNCT(FILNAM) 1363c WRITE(6,*) 'open file' 1364c: WRITE(6,*) filnam 1365 call ccpdpn(3,FILNAM,'READONLY','F',0,0) 1366C 1367 NATM = 1 13685 READ(3,10,ERR=5,END=20) HEAD, ATOM(NATM), SEQ(NATM), RES1 1369 1 , (X(I,NATM),I=1,3),occ(natm),bfac(natm) 1370c normal format 1371c10 FORMAT(A6,7X,2A4,A5,4X,3F8.3,2f6.2) 1372c xplor format 137310 FORMAT(A6,6X,2A4,1X,A5,4X,3F8.3,2f6.2) 1374 IF (HEAD.NE.'ATOM '.AND.HEAD.NE.'HETATM') GOTO 5 1375 IF (ATOM(NATM)(1:1).EQ.'H'.or.ATOM(NATM)(1:2).EQ.' H') goto 5 1376C 1377 IF (RES1(1:1).NE.' '.AND.RES1(2:2).EQ.' ') THEN 1378 RES2(1:1) = RES1(1:1) 1379 RES2(2:4) = RES1(3:5) 1380 ELSE 1381 RES2 = RES1(2:5) 1382 END IF 1383C 1384 call spstrunct(atom(natm)) 1385 call spstrunct(seq(natm)) 1386 call spstrunct(RES2) 1387 res(natm)=res2 1388C 138912 CONTINUE 1390C 1391 NATM = NATM + 1 1392 IF (NATM.GT.MAXATM) STOP 'Atoms can not be more than MAXATM.' 1393 GOTO 5 139420 CONTINUE 1395 NATM = NATM - 1 1396 CLOSE (3) 1397C 1398 END 1399c 1400c 1401 SUBROUTINE GETPDB2(X,ATOM,SEQ,RES,BFAC,OCC,NATM,FILNAM) 1402C A subroutine to read the Protein data bank format coordinate file 1403c X 3*maxatom arrays, output xyz coordinates 1404c ATOM maxatom character*4 arrays, output atom names 1405c RES maxatom character*4 arrays, output residue number e.g. A101 A103... 1406C If residue number is like A 1 then it will become A1 1407c SEQ maxatom character*4 arrays, output peptide name. e.g. PHE TRY .... 1408C NATM is the number of atoms in the coordinate file 1409c FILNAT is the file name of the coordinate. 1410c 1411c in this subroutine, hydrogen atoms were included. 1412c 1413c getpdb1 is for xplor file and exclude Hydrogen 1414c getpdb2 is for any file 1415 parameter (maxatom = 50000) 1416 DIMENSION X(3,maxatom) 1417 real bfac(maxatom),occ(maxatom) 1418 CHARACTER*4 ATOM(maxatom),RES(maxatom),SEQ(maxatom),RES2 1419 CHARACTER FILNAM*80,HEAD*6,RES1*5 1420C DATA FILNAM/'DIAMOND'/ 1421C 1422 CALL SPSTRUNCT(FILNAM) 1423c WRITE(6,*) 'open file' 1424c: WRITE(6,*) filnam 1425 call ccpdpn(3,FILNAM,'READONLY','F',0,0) 1426C 1427 NATM = 1 14285 READ(3,10,ERR=5,END=20) HEAD, ATOM(NATM), SEQ(NATM), RES1 1429 1 , (X(I,NATM),I=1,3),occ(natm),bfac(natm) 1430c normal format 1431c10 FORMAT(A6,7X,2A4,A5,4X,3F8.3,2f6.2) 1432c xplor format 143310 FORMAT(A6,6X,2A4,1X,A5,4X,3F8.3,2f6.2) 1434 IF (HEAD.NE.'ATOM '.AND.HEAD.NE.'HETATM') GOTO 5 1435c IF (ATOM(NATM)(1:1).EQ.'H'.or.ATOM(NATM)(1:2).EQ.' H') goto 5 1436C 1437 IF (RES1(1:1).NE.' '.AND.RES1(2:2).EQ.' ') THEN 1438 RES2(1:1) = RES1(1:1) 1439 RES2(2:4) = RES1(3:5) 1440 ELSE 1441 RES2 = RES1(2:5) 1442 END IF 1443C 1444 call spstrunct(atom(natm)) 1445 call spstrunct(seq(natm)) 1446 call spstrunct(RES2) 1447 res(natm)=res2 1448C 144912 CONTINUE 1450C 1451 NATM = NATM + 1 1452 IF (NATM.GT.maxatom) STOP 'Atoms can not be more than maxatom.' 1453 GOTO 5 145420 CONTINUE 1455 NATM = NATM - 1 1456 CLOSE (3) 1457C 1458 END 1459c 1460c 1461 SUBROUTINE GETPDB3(X,ATOM,SEQ,CH,RES,BFAC,OCC,NATM,FILNAM) 1462C A subroutine to read the Protein data bank format coordinate file 1463c X 3*15000 arrays, output xyz coordinates 1464c ATOM 15000 character*4 arrays, output atom names 1465c RES 15000 character*4 arrays, output residue number e.g. A101 A103... 1466C If residue number is like A 1 then it will become A1 1467c SEQ 15000 character*4 arrays, output peptide name. e.g. PHE TRY .... 1468C NATM is the number of atoms in the coordinate file 1469c FILNAT is the file name of the coordinate. 1470c 1471c in this subroutine, hydrogen atoms were included. 1472c 1473c getpdb1 is for xplor file and exclude Hydrogen 1474c getpdb2 is for any file 1475 parameter (maxpro = 50000) 1476 DIMENSION X(3,maxpro) 1477 real bfac(maxpro),occ(maxpro) 1478 CHARACTER*4 ATOM(maxpro),RES(maxpro),SEQ(maxpro),RES2 1479 character*1 ch(maxpro) 1480 CHARACTER FILNAM*80,HEAD*6,RES1*5 1481C 1482 CALL SPSTRUNCT(FILNAM) 1483c WRITE(6,*) 'open file' 1484c: WRITE(6,*) filnam 1485 call ccpdpn(3,FILNAM,'READONLY','F',0,0) 1486C 1487 NATM = 1 14885 READ(3,10,ERR=5,END=20) HEAD, ATOM(NATM), SEQ(NATM), ch(natm), 1489 1 res(natm), (X(I,NATM),I=1,3),occ(natm),bfac(natm) 1490c normal format 1491c10 FORMAT(A6,7X,2A4,A5,4X,3F8.3,2f6.2) 1492c xplor format 149310 FORMAT(A6,6X,2A4,1X,a1,A4,4X,3F8.3,2f6.2) 1494 IF (HEAD.NE.'ATOM '.AND.HEAD.NE.'HETATM') GOTO 5 1495c IF (ATOM(NATM)(1:1).EQ.'H'.or.ATOM(NATM)(1:2).EQ.' H') goto 5 1496C 1497c IF (RES1(1:1).NE.' '.AND.RES1(2:2).EQ.' ') THEN 1498c RES2(1:1) = RES1(1:1) 1499c RES2(2:4) = RES1(3:5) 1500c ELSE 1501c RES2 = RES1(2:5) 1502c END IF 1503C 1504 call spstrunct(atom(natm)) 1505 call spstrunct(seq(natm)) 1506c call spstrunct(RES2) 1507c res(natm)=res2 1508C 150912 CONTINUE 1510C 1511 NATM = NATM + 1 1512 IF (NATM.GT.maxpro) then 1513 WRITE(6,*) 'Error: atom number is',natm 1514 WRITE(6,*) 'Error.. Atoms cannot be more than',maxpro 1515 stop 1516 end if 1517 GOTO 5 151820 CONTINUE 1519 NATM = NATM - 1 1520 CLOSE (3) 1521C 1522 END 1523c 1524c 1525 SUBROUTINE GETPDBCA(X,ATOM,RES,SEQ,NATM,FILNAM) 1526C A subroutine to read the Protein data bank format coordinate file 1527c In this version only Ca atoms are read. -- Guoguang 950922 1528c X 3*50000 arrays, output xyz coordinates 1529c ATOM 50000 character*4 arrays, output atom names 1530c RES 50000 character*4 arrays, output residue number e.g. A101 A103... 1531C If residue number is like A 1 then it will become A1 1532c SEQ 50000 character*4 arrays, output peptide name. e.g. PHE TRY .... 1533C NATM is the number of atoms in the coordinate file 1534c FILNAT is the file name of the coordinate. 1535c 1536 DIMENSION X(3,50000) 1537 CHARACTER*4 ATOM(50000),RES(50000),SEQ(50000),RES2 1538 CHARACTER FILNAM*80,HEAD*6,RES1*5 1539C DATA FILNAM/'DIAMOND'/ 1540C 1541 CALL SPSTRUNCT(FILNAM) 1542 call ccpdpn(3,FILNAM,'READONLY','F',0,0) 1543C 1544 NATM = 1 15455 READ(3,10,ERR=5,END=20) HEAD, ATOM(NATM), SEQ(NATM), RES1 1546 1 , (X(I,NATM),I=1,3) 154710 FORMAT(A6,7X,2A4,A5,4X,3F8.3) 1548 IF (HEAD(1:3).EQ.'END') GOTO 20 1549 IF (HEAD.NE.'ATOM '.AND.HEAD.NE.'HEDATM') GOTO 5 1550 IF (ATOM(NATM).NE.'CA ') GOTO 5 1551C 1552 IF (RES1(1:1).NE.' '.AND.RES1(2:2).EQ.' ') THEN 1553 RES2(1:1) = RES1(1:1) 1554 RES2(2:4) = RES1(3:5) 1555 ELSE 1556 RES2 = RES1(2:5) 1557 END IF 1558C 1559 call spstrunct(atom(natm)) 1560 call spstrunct(seq(natm)) 1561 call spstrunct(RES2) 1562 res(natm)=res2 1563C 156412 CONTINUE 1565C 1566 NATM = NATM + 1 1567 IF (NATM.GT.50000) STOP 'Atoms cannot be more than 50000.' 1568 GOTO 5 156920 CONTINUE 1570 NATM = NATM - 1 1571 CLOSE (3) 1572C 1573 END 1574c 1575c 1576 SUBROUTINE GETRDI(X,ATOM,RES,SEQ,NATM,FILNAM) 1577C A subroutine to read the Diamond format coordinate file 1578c X 3*5500 arrays, output xyz coordinates 1579c ATOM 5500 character*4 arrays, output atom names 1580c RES 5500 character*4 arrays, output residue number e.g. A101 A103... 1581C If res is 'A 1' then it will be truncated to 'A1 ' 1582c SEQ 5500 character*4 arrays, output peptide name. e.g. PHE TRY .... 1583C NATM is the number of atoms in the coordinate file 1584c FILNAT is the file name of the coordinate. 1585c 1586 DIMENSION X(3,5500) 1587 CHARACTER*4 ATOM(5500),RES(5500),SEQ(5500),RES1*4 1588 CHARACTER FILNAM*80 1589C DATA FILNAM/'DIAMOND'/ 1590C 1591 CALL SPSTRUNCT(FILNAM) 1592 call ccpdpn(2,FILNAM,'READONLY','F',0,0) 1593C 1594 READ(2,'(A)') ATOM(1),ATOM(1),ATOM(1) 1595C 1596 NATM = 1 15975 READ(2,10,ERR=5,END=20) (X(I,NATM),I=1,3), 1598 1 SEQ(NATM), RES(NATM), ATOM(NATM) 159910 FORMAT(3F10.4,35X,A3,A4,3X,A4) 1600 IF (RES(NATM).EQ.'END '.OR.RES(natm).EQ.'CAST') GOTO 5 1601C 1602 CALL SPSTRUNCT(ATOM(NATM)) 1603 CALL SPSTRUNCT(SEQ(NATM)) 1604 CALL SPSTRUNCT(RES(NATM)) 1605C 160612 NATM = NATM + 1 1607 GOTO 5 160820 CONTINUE 1609 NATM = NATM - 1 1610 I = 1 1611 CLOSE (2) 1612C 161325 IF (I.LE.NATM.AND.SEQ(I).NE.' ') THEN 1614C 1615 DO J = I-1, 1, -1 1616 IF (J.LT.1) GOTO 30 1617 IF (RES(I).EQ.RES(J)) THEN 1618 SEQ(J) = SEQ(I) 1619 ELSE 1620 GOTO 30 1621 END IF 1622 END DO 1623C 162430 DO J = I+1, NATM 1625 IF (RES(I).EQ.RES(J)) THEN 1626 SEQ(J) = SEQ(I) 1627 ELSE 1628 GOTO 40 1629 END IF 1630 END DO 1631C 163240 CONTINUE 1633 I = J 1634 GOTO 25 1635 ELSE IF (I.LT.NATM) THEN 1636 I = I + 1 1637 GOTO 25 1638 ELSE 1639 RETURN 1640 END IF 1641 1642 END 1643c 1644c 1645 SUBROUTINE GETSHELX ( X, ATOM, NATM, FILNAM ) 1646c A routine to get coordinates and atom names from a SHELXX format 1647c file. 1648c X is 3*500 array to save the coordinates. 1649c ATOM is 4 bytes 500-dim array to represent the atom names 1650c NATM represents the number of atoms 1651c FILNAM is the SHELXX format file name. 1652c Sample coordinate format: 1653c FVAR 0.90195 1654c ND 5 0.59828 0.24644 10.24980 11.00000 0.01625 0.00675 = 1655c 0.03496 0.00372 0.00326 0.00215 1656c O1 4 0.58560 0.24368 0.55481 11.00000 0.02699 0.04010 = 1657c 0.01761 0.00117 0.00328 0.01509 1658c ...... 1659c ...... 1660c END 1661c 1662 PARAMETER (MAXATM = 500 ) 1663 CHARACTER FILNAM*80,ATOM(MAXATM)*4,KEY*80 1664 DIMENSION X(3,MAXATM) 1665C 1666C Open the SHELXX format coordinate file. 1667C 1668 call ccpdpn(1,FILNAM,'READONLY','F',0,0) 1669C 1670 NATM = 0 1671 READ(1,*) 1672C 16735 NATM = NATM + 1 1674 IF (NATM.GT.MAXATM) 1675 + STOP 'Atoms of Shelxx cannot be more than 500.' 1676 READ(1,10) ATOM(NATM),(X(I,NATM),I=1,3) 1677C WRITE(6,*) NATM,ATOM(NATM),(X(I,NATM),I=1,3) 167810 FORMAT(A4,5X,3F10.5) 1679 IF (ATOM(NATM).EQ.'END '.OR. ATOM(NATM).EQ.' ') GOTO 20 1680 READ(1,*) 1681 GOTO 5 168220 NATM = NATM -1 1683 CLOSE (1) 1684 RETURN 1685 END 1686C 1687C * ROH - MATRIX 1688C 1689c 1690 SUBROUTINE HUBER(ROT,E) 1691 1692 DIMENSION E(3,3),ROT(3) 1693c EQUIVALENCE (PS,AL,PA), (TH,BE,PB), (PH,GA,PC) 1694c EQUIVALENCE (S1,SIPS), (S2,SITH), (S3,SIPH) 1695c 1 ,(C1,COPS), (C2,COTH), (C3,COPH) 1696 1697 PS=ROT(1) 1698 TH=ROT(2) 1699 PH=ROT(3) 1700 1701 ARC=3.14159265359/180.0 1702 SIPS = SIN(PS*ARC) 1703 COPS = COS(PS*ARC) 1704 SITH = SIN(TH*ARC) 1705 COTH = COS(TH*ARC) 1706 SIPH = SIN(PH*ARC) 1707 COPH = COS(PH*ARC) 1708 E(1,1) = COPS*COPH - SIPS*SIPH*SITH 1709 E(1,2) =-SIPS*COTH 1710 E(1,3) = COPS*SIPH + SIPS*SITH*COPH 1711 E(2,1) = SIPS*COPH + COPS*SITH*SIPH 1712 E(2,2) = COPS*COTH 1713 E(2,3) = SIPS*SIPH - COPS*SITH*COPH 1714 E(3,1) =-COTH*SIPH 1715 E(3,2) = SITH 1716 E(3,3) = COTH*COPH 1717 END 1718C 1719C 1720C * ROH - MATRIX 1721C 1722 SUBROUTINE HUBERARC(ROT,E) 1723 1724 DIMENSION E(3,3),ROT(3) 1725 EQUIVALENCE (PS,AL,PA), (TH,BE,PB), (PH,GA,PC) 1726 EQUIVALENCE (S1,SIPS), (S2,SITH), (S3,SIPH) 1727 1 ,(C1,COPS), (C2,COTH), (C3,COPH) 1728 1729 PS=ROT(1) 1730 TH=ROT(2) 1731 PH=ROT(3) 1732 1733C ARC=3.14159265359/180.0 1734C ARC= 1.0 1735 SIPS = SIN(PS) 1736 COPS = COS(PS) 1737 SITH = SIN(TH) 1738 COTH = COS(TH) 1739 SIPH = SIN(PH) 1740 COPH = COS(PH) 1741 E(1,1) = COPS*COPH - SIPS*SIPH*SITH 1742 E(1,2) =-SIPS*COTH 1743 E(1,3) = COPS*SIPH + SIPS*SITH*COPH 1744 E(2,1) = SIPS*COPH + COPS*SITH*SIPH 1745 E(2,2) = COPS*COTH 1746 E(2,3) = SIPS*SIPH - COPS*SITH*COPH 1747 E(3,1) =-COTH*SIPH 1748 E(3,2) = SITH 1749 E(3,3) = COTH*COPH 1750 END 1751c 1752c 1753 SUBROUTINE IARRGIVE(N,XYZ0,XYZ) 1754C Copy an N-dimensional integer array XYZ0 to another one XYZ. 1755 INTEGER*4 XYZ(N),XYZ0(N) 1756 DO I = 1, N 1757 XYZ(I) = XYZ0(I) 1758 END DO 1759 END 1760c 1761c 1762C------MULTIPLY AN INTEGER ARRAY BY A CONSTANT 1763 SUBROUTINE IARRMC(IM,IN,A1,C,A) 1764 INTEGER A1(IM,IN),A(IM,IN),C 1765 DO 10 I1=1,IM 1766 DO 10 I2=1,IN 176710 A(I1,I2)=A1(I1,I2)*C 1768 END 1769c 1770c 1771 subroutine imatext(isym0,sout,ss1,ss2,ss3,trans) 1772c A subroutine to transfer a symmetric matrix into a string. 1773c isym0 is 3*3 symmetric matrix. 1774c in is 3-dimensional strings of the input like 'x''y''z' or 'h''k''l' 1775c trans is true. recognize the matrix as a transpose one. 1776c sout is a 3-dimensional output string 1777c 1778c Guoguang 91-02-06 1779 character sout(3)*12,ss1*1,ss2*1,ss3*1 1780 logical*1 trans 1781 integer isym0(3,3) 1782 integer isym(3,3) 1783 character out(3,3)*4,temp*4,in(3)*1 1784c equivalence (sout,out) 1785c 1786 in(1) = ss1 1787 in(2) = ss2 1788 in(3) = ss3 1789C 1790 if (trans) then 1791 do j = 1,3 1792 do i = 1,3 1793 isym(i,j) = isym0(j,i) 1794 end do 1795 end do 1796 else 1797 do j = 1, 3 1798 do i = 1, 3 1799 isym(i,j) = isym0(i,j) 1800 end do 1801 end do 1802 end if 1803c 1804 do i = 1,3 1805 do j = 1,3 1806 out(i,j)(4:4) = in(j) 1807 if (isym(i,j).eq.0) then 1808 out(i,j)(1:4) = ' ' 1809 else if(abs(isym(i,j)).ge.10) then 1810 write(6,*) 'isym(',i,j,')=',isym(i,j),' > 10' 1811 stop 'error' 1812 else 1813 if (isym(i,j).gt.0) out(i,j)(1:1) = '+' 1814 if (isym(i,j).lt.0) out(i,j)(1:1) = '-' 1815 if (abs(isym(i,j)).ne.1) then 1816 write(out(i,j)(2:3),'(i1,a1)') abs(isym(i,j)),'*' 1817 else 1818 out(i,j)(2:3) = ' ' 1819 end if 1820 end if 1821 end do 1822 end do 1823c 1824 do i = 1, 3 1825 write(sout(i),'(3a4)') (out(i,j),j=1,3) 1826c write(6,*) sout(i) 1827c if (sout(i)(1:1).eq.'+') sout(i)(1:1) = ' ' 1828c make the text shorter 1829 il = lenstr(sout(i)) 1830 if (il.gt.1) then 1831 iln = il 1832 do j = il-1, 1, -1 1833 if (sout(i)(j:j).eq.' ') then 1834 lg = iln - j 1835 sout(i)(j:iln-1) = sout(i)(j+1:iln) 1836 sout(i)(iln:iln) = ' ' 1837 iln = iln - 1 1838 end if 1839 end do 1840 end if 1841 if (sout(i)(1:1).eq.'+') sout(i)(1:1) = ' ' 1842 end do 1843c 1844 end 1845c 1846c 1847C-------A STANDARD SUBPROGRAM TO MULTIPLY TWO ARRAYS (a new quick version 1848C --- this is for I4 version 1849C) BMC 91/02/11/ Guoguang 1850C A1(IM1,IN1)*A2(IM2,IN2)=RSLT(IM1,IN2) 1851 SUBROUTINE IMATMULT(IM1,IN1,IM2,IN2,A1,A2,RSLT) 1852 INTEGER A1(IM1,IN1),A2(IM2,IN2),RSLT(IM1,IN2) 1853 IF (IN1.NE.IM2) THEN 1854 stop 'The two arrays cannot be multiplied' 1855 end if 1856 DO I = 1, IM1 1857 DO J = 1, IN2 1858 RSLT(I,J) = 0. 1859 DO k = 1, in1 1860 RSLT(I,J) = RSLT(I,J) + A1(I,K)*A2(K,J) 1861 END DO 1862 END DO 1863 END DO 1864 RETURN 1865 END 1866c 1867c 1868c if you have a rigid movement x2=R*x1+1 1869c where R(3,3) is the rotation and t(3) is the translation 1870c then x1 =R2*x1-t2 1871c where R2=R-1 1872c t2 =R-1*t 1873c The program reads in R,T matrix (3,4) and returns a inversed RT2 1874c nmat is the number of matrix 1875c 1876 subroutine invrt(nmat,rt,rt2) 1877 real rt(3,4,nmat),rt2(3,4,nmat) 1878 real b1(3),b2(3) 1879 integer me(3) 1880c 1881 do i = 1, nmat 1882 call arrgive(9,rt(1,1,i),rt2(1,1,i)) 1883 call ivsn(3,rt2(1,1,i),b1,b2,me,de,1.e-4) 1884 call arrmc(3,1,rt(1,4,i),-1.,b1) 1885 call matmult(3,3,3,1,rt2(1,1,i),b1,rt2(1,4,i)) 1886 end do 1887c 1888 end 1889c 1890c 1891C--------- IN PLACE MATRIX INVERSION --------- 1892 1893 SUBROUTINE IVSN(N,A,B,C,ME,DE,EP) 1894 1895C---- N is the dimension of the matrix A with N rows and N columns 1896c---- A-N,N- is the input matrix with N rows and N columns. 1897c--- ME(N), B(N) AND C(N) are the buffer of N-dimensional array. 1898c--- DE could any value 1899c EP is the least value of the pivot i.e. when |A| < EP, the program 1900c will stop. 1901 1902c on output, DE is the determinant of the original matrix 1903c on output, A is the inverse of the original input A 1904c ME is an index array which keeps track of row swaps 1905 1906 DIMENSION A(N,N),ME(N),B(N),C(N) 1907 INTEGER I,J,K,ME 1908 1909 DE = 1.0 1910 1911 DO 10 J=1,N 191210 ME(J)=J 1913 1914 DO 20 I=1,N 1915 Y=0. 1916 DO 30 J=I,N 1917 IF (ABS(A(I,J)).LE.ABS(Y)) GOTO 30 1918 K=J 1919 Y=A(I,J) 192030 CONTINUE 1921 DE = DE*Y 1922 IF (ABS(Y).LT.EP) then 1923 write(6,*) 'Ill conditioned matrix:' 1924 write(6,'(3f12.6)') a 1925 WRITE(6,*) 'the pivot of the matrix is ', y, ' N = ',N 1926 STOP 4444 1927 end if 1928 Y = 1.0/Y 1929 DO 40 J=1,N 1930 C(J) = A(J,K) 1931 A(J,K) = A(J,I) 1932 A(J,I) = -C(J)*Y 1933 B(J) = A(I,J)*Y 193440 A(I,J) = A(I,J)*Y 1935 A(I,I) = Y 1936 J = ME(I) 1937 ME(I) = ME(K) 1938 ME(K) = J 1939 DO 11 K=1,N 1940 IF(K.EQ.I) GOTO 11 1941 DO 12 J=1,N 1942 IF (J.EQ.I) GOTO 12 1943 A(K,J) = A(K,J) - B(J)*C(K) 194412 CONTINUE 194511 CONTINUE 1946 194720 CONTINUE 1948 1949 DO 33 I=1,N 1950 DO 44 K=1,N 1951 IF(ME(K).EQ.I) GOTO 55 195244 CONTINUE 195355 IF(K.EQ.I) GOTO 33 1954 DO 66 J=1,N 1955 W = A(I,J) ! swap rows k and i 1956 A(I,J) = A(K,J) 195766 A(K,J) = W 1958 IW = ME(I) ! keep track of index changes 1959 ME(I) = ME(K) 1960 ME(K) = IW 1961 DE = -DE ! determinant changes sign when two rows swapped 196233 CONTINUE 1963 RETURN 1964 END 1965c 1966c 1967 SUBROUTINE KABMOD(TH1,TH2,TH3,TH11,TH12,TH13) 1968 DIMENSION TH(6) 1969C 1970C --- ERRECHNET EULERWINKEL IM BEREICH -180<TH<180 1971C 1972 TH(1) = TH1 1973 TH(2) = TH2 1974 TH(3) = TH3 1975 TH(4) = TH11 1976 TH(5) = TH12 1977 TH(6) = TH13 1978 DO 1 I=1,6 1979 IF (TH(I).LE.-180.) TH(I) = TH(I) + 360. 1980 1 IF (TH(I).GT.180.) TH(I) = TH(I) - 360. 1981 TH1 = TH(1) 1982 TH2 = TH(2) 1983 TH3 = TH(3) 1984 TH11 = TH(4) 1985 TH12 = TH(5) 1986 TH13 = TH(6) 1987 RETURN 1988 END 1989c 1990c 1991 SUBROUTINE LATTIC(NSYM,SYM,LAT) 1992C A subroutine to add the translation vector to the symmetry arrays 1993c and modify the number of symmetry operations according to the 1994c Bravais Lattice. 1995C CHARACTER*1 LAT 1996C DIMENSION SYM(3,4,192) 1997c Parameter description: 1998c NSYM are the number of the symmetric operation. The input value do 1999c not include the non-1 translation operation. But output includes it. 2000c SYM(3,4,NSYM) is the symmetry op arrays. 2001c LAT is the Bravais lattice. It could be "P", "A", "B", "C", "I", 2002c "F",or "R" 2003C Written by Guoguang Lu 2004c 27/02/88 2005 CHARACTER*1 LAT 2006 DIMENSION SYM(3,4,192),TR(3,4) 2007c 2008 IF (LAT.EQ.'P') RETURN 2009 CALL ARRMC(3,4,TR,0.,TR) 2010 IF (LAT.EQ.'A') THEN 2011 NTR = 2 2012 TR(2,2) = .5 2013 TR(3,2) = .5 2014 ELSE IF (LAT.EQ.'B') THEN 2015 NTR = 2 2016 TR(1,2) = .5 2017 TR(3,2) = .5 2018 ELSE IF (LAT.EQ.'C') THEN 2019 NTR = 2 2020 TR(2,2) = .5 2021 TR(1,2) = .5 2022 ELSE IF (LAT.EQ.'I') THEN 2023 NTR = 2 2024 TR(1,2) = .5 2025 TR(2,2) = .5 2026 TR(3,2) = .5 2027 ELSE IF (LAT.EQ.'F') THEN 2028 NTR = 4 2029 TR(1,2) = .5 2030 TR(2,2) = .5 2031 TR(2,3) = .5 2032 TR(3,3) = .5 2033 TR(3,4) = .5 2034 TR(1,4) = .5 2035 ELSE IF (LAT.EQ.'R') THEN 2036 NTR = 3 2037 TR(1,2) = 1./3. 2038 TR(2,2) = 2./3. 2039 TR(3,2) = 2./3. 2040 TR(1,3) = 2./3. 2041 TR(2,3) = 1./3. 2042 TR(3,3) = 1./3. 2043 ELSE 2044 WRITE(6,*) 'Error Lattice>> No lattice ',LAT 2045 WRITE(6,*) 'Only P,A,B,C,I,F and R are allowed.' 2046 END IF 2047 DO ISYM = 1, NSYM 2048 DO ITR = 2, NTR 2049 I = (ITR-1)*NSYM + ISYM 2050 CALL ARRMC(3,3, SYM(1,1,ISYM), 1., SYM(1,1,I) ) 2051 CALL ARRAD(3,1, TR(1,ITR), SYM(1,4,ISYM), SYM(1,4,I) ) 2052 DO J = 1, 3 20535 IF (SYM(J,4,I).GE.1.) THEN 2054 SYM(J,4,I) = SYM(J,4,I) - 1. 2055 GOTO 5 2056 END IF 20576 IF (SYM(J,4,I).LT.-1.) THEN 2058 SYM(J,4,I) = SYM(J,4,I) + 1. 2059 GOTO 6 2060 END IF 2061 END DO ! J 2062 END DO ! ITR 2063 END DO ! ISYM 2064 NSYM = NSYM * NTR 2065 RETURN 2066 END 2067c 2068c 2069 SUBROUTINE LSQEQ(M,N,A,PHI,DETA,ATA,B1) 2070C DIMENSION A(M,N),PHI(M),DETA(N),ATA(N,N) 2071C DIMENSION B1(N) 2072C This is a subroutine to compute the solution of least square equation 2073C of a normal matrix. i.e. AT*A*DETA = AT*PHI ==> DETA = (AT*A)^(-1)*AT*PHI 2074C It could be get the deta the make the SIGMA (phi)^2 = minimum 2075c where there are M equations 2076c ....... 2077c phi i = 0 2078c ....... 2079c 2080c Note: In the program parameter PHI(I) = - phi0 i in the equation. 2081c 2082C written by Guoguang Lu 2083c 23/01/1989 2084c parameter description: 2085c M is the number of equations. i.e. PHIi = 0 (i=1,...M) 2086C N is the number of the variables. 2087C (N should be less than 3*N0 given in parameter statement.) 2088C (At this moment, N0 is 6000.) 2089c (A,PHI) is a M*(N+1) normal matrix 2090c A is a M*N array 2091c PHI is a M-dimensional vector. Here, PHI should be -phi0 in the equation. 2092c DETA is the output shift. 2093c ATA(N,N ) is a N*N buffer array. 2094c B1(N) is a N-dimensional buffer array. 2095c Note : M must not be less than N. 2096 PARAMETER (N0 = 6000 ) 2097 DIMENSION A(M,N),PHI(M),DETA(N),ATA(N,N) 2098 DIMENSION B1(N),ME(3*N0) 2099 IF (M.LT.N) THEN 2100 STOP 'Equation number is not enough' 2101 else if (M.EQ.N) THEN 2102 CALL IVSN(M,A,B1,DETA,ME,DE,1.E-5) 2103c CALL ARRMC(N,1,DETA,0.,DETA) 2104 call arrvalue(n,deta,0.) 2105 DO 50 I = 1, N 2106 DO 50 J = 1, N 210750 DETA(I) = DETA(I) + ATA(I,J)*PHI(J) 2108 ELSE 2109C 2110C T 2111C ATA = A * A 2112C 2113 call arrvalue(n*n,ata,0.) 2114c CALL ARRMC(N,N,ATA,0.,ATA) 2115 DO 10 I = 1,N 2116 DO 10 J = I,N 2117 DO 10 K = 1,M 211810 ATA(I,J) = ATA(I,J) + A(K,I)*A(K,J) 2119 IF (N.GE.2) THEN 2120 DO 20 I = 1, N 2121 DO 20 J = 1, I-1 212220 ATA(I,J) = ATA(J,I) 2123 END IF 2124C 2125C 2126C ATA = ATA^-1 2127C 2128 CALL IVSN (N,ATA,B1,DETA,ME,DE,1.E-5) 2129C 2130C T 2131C B1 = A * PHI 2132C 2133c CALL ARRMC(N,1,B1,0.,B1) 2134 call arrvalue(n,b1,0.) 2135 DO 30 I = 1, N 213630 B1(I) = POIMULT(M, M, A(1,I), PHI) 2137C 2138C DETA = ATA^-1 * B1 2139C 2140C T T 2141C i.e. DETA = (A * A) * A * PHI 2142C 2143c CALL ARRMC(N,1,DETA,0.,DETA) 2144 call arrvalue(n,deta,0.) 2145 DO 40 I = 1, N 2146 DO 40 J = 1, N 214740 DETA(I) = DETA(I) + ATA(I,J)*B1(J) 2148C 2149 END IF 2150C 2151 END 2152c 2153c 2154 subroutine matitor(im,in,imat,rmat) 2155c Convert a matrix from integer*2 to real*4 2156c 2157 real*4 rmat(im,in) 2158 integer imat(im,in) 2159c 2160 do j = 1, in 2161 do i = 1, im 2162 rmat(i,j) = imat(i,j) 2163 end do 2164 end do 2165 end 2166c 2167c 2168C-------A STANDARD SUBPROGRAM TO MULTIPLY TWO ARRAYS (a new quick version 2169C) BMC 89/11/05/ Guoguang 2170C A1(IM1,IN1)*A2(IM2,IN2)=RSLT(IM1,IN2) 2171 SUBROUTINE MATMULT(IM1,IN1,IM2,IN2,A1,A2,RSLT) 2172 DIMENSION A1(IM1,IN1),A2(IM2,IN2),RSLT(IM1,IN2) 2173 IF (IN1.NE.IM2) THEN 2174 stop 'The two arrays cannot be multiplied' 2175 end if 2176 DO I = 1, IM1 2177 DO J = 1, IN2 2178 RSLT(I,J) = 0. 2179 DO k = 1, in1 2180 RSLT(I,J) = RSLT(I,J) + A1(I,K)*A2(K,J) 2181 END DO 2182 END DO 2183 END DO 2184 RETURN 2185 END 2186c 2187c 2188 subroutine matrtoi(im,in,rmat,imat) 2189c Convert a matrix from real to integer 2190c 2191 real*4 rmat(im,in) 2192 integer imat(im,in) 2193c 2194 do j = 1, in 2195 do i = 1, im 2196 imat(i,j) = rmat(i,j) 2197 end do 2198 end do 2199 end 2200c 2201c 2202 FUNCTION MATSYM(S,TEXT,ICOL) 2203C 2204C A FORTRAN FUNCTION SUBPROGRAM TO SCAN A SYMMETRY CARD AND BUILD 2205C THE CORRESPONDING SYMMETRY-OPERATION MATRIX. THE CONTENTS OF THE 2206C CARD AND ERRORS ENCOUNTERED ARE WRITTEN OFF-LINE IN THE CALLING 2207C PROGRAM. THE FUNCTIONAL VALUE IS NORMALLY ZERO, WHEN AN ERROR HAS 2208C BEEN DETECTED THE NUMBERS 1-4 ARE GIVEN BACK: 2209C ERROR-NUMBER 1: INVALID CHARACTER AT POSITION .... 2210C 2: SYNTAX ERROR AT NONBLANK POSITION .... 2211C 3: BAD COMMAS 2212C 4: DETERMINANT OF ROTATION MATRIX NOT +1 OR -1 2213C 2214C EXPLANATION OF ARGUMENTS 2215C S IS THE 1ST ELEMENT OF THE 3X4 ARRAY WHERE THE SYM. OP. IS TO BE 2216C STORED. 2217C TEXT IS THE ARRAY OF THE TEXT RECORD WHICH IS TO BE SCANNED (LENGTH: 2218C 36 CHARACTERS, CONTIGUOUSLY) 2219C ICOL POINTS - IN CASE OF ERRORS - TO THE POSITION OF AN INVALID 2220C CHARACTER OR AN SYNTAX ERROR WITHIN THE SYMMETRY CARD 2221C 2222C MODIFIED FOR PDP-11 06/JAN/78 S.J. REMINGTON 2223C ***** LAST UPDATE 07/10/76 MRS WRS:PR2.MATSYM 2224C ***** PREVIOUS UPDATES 07/07/75 06/11/72 2225C 2226C 2227 DIMENSION S(3,4), ICHAR(36) 2228 CHARACTER*1 CHAR,VALID(14),BLANK 2229 CHARACTER*36 TEXT 2230 EQUIVALENCE (BLANK,VALID(14)) 2231 DATA VALID/'1','2','3','4','5','6','X','Y','Z','-','+','/',',', 2232 . ' '/ 2233C 2234C----INITIALIZE SIGNALS 2235 ICOL = 0 2236 IFIELD = 0 2237C----CLEAR SYMOP ARRAY 2238 DO 4 J=1,4 2239 DO 4 I=1,3 2240 4 S(I,J) = 0.0 2241C 2242C----TEST THE SYMMETRY TEXT (TEXT) FOR ILLEGAL CHARACTERS 2243C 2244C THE STATEMENT IS SCANNED AND ARRAY ICHAR IS FILLED WITH THE INDICES 2245C OF EACH ELEMENT OF TEXT INTO THE ARRAY 'VALID' (I.E. IF THE SEVENTH 2246C CHAR OF 'TEXT' -IGNORING BLANKS- IS 'Y', ICHAR(7)=8) 2247C 2248 ICOLMX = 0 2249 DO 20 I=1,36 2250 CHAR = TEXT(I:I) 2251 IF(CHAR .EQ. BLANK) GO TO 20 2252 ICOLMX = ICOLMX + 1 2253 DO 10 J=1,13 2254 IF( VALID(J) .NE. CHAR) GO TO 10 2255 ICHAR(ICOLMX) = J 2256 GO TO 20 2257 10 CONTINUE 2258C--IF GOT TO HERE, THE CHAR ISN'T IN 'VALID' 2259 JTXTSC = I 2260 GO TO 9975 2261 20 CONTINUE 2262C 2263C----BEGIN FIELD LOOP 2264C 2265 101 IFIELD = IFIELD + 1 2266 IF(IFIELD-3) 104,104,103 2267C HAVE ALL CHARACTERS BEEN ANALYSED IN 3 FIELDS? 2268 103 IF(ICOL-ICOLMX) 9987,1000,1000 2269C---NO MORE THAN THREE FIELDS PERMITTED 2270 104 SIGN = 1.0 2271 ICOL = ICOL+1 2272C----MARCH ON, FIND WHAT'S NEXT 2273 IFNUM=ICHAR(ICOL) 2274 IF (11-IFNUM) 9980,110,110 2275C----TEST FOR SIGNS 2276 110 IF (10-IFNUM) 118,115,112 2277C----TEST TO DISTINGUISH BETWEEN NUMBERS AND LETTERS 2278 112 IF (6-IFNUM) 125,130,130 2279C----FOUND A SIGN 2280C COME HERE FOR MINUS 2281 115 SIGN = -1.0 2282C COME HERE FOR PLUS 2283 118 ICOL = ICOL + 1 2284C----NEXT CHARACTER MUST BE A NUMBER OR LETTER 2285 IFNUM=ICHAR(ICOL) 2286 IF (IFNUM-9) 122,122,9980 2287 122 IF (IFNUM-6) 130,130,125 2288C----A LETTER 2289 125 IFNUM = IFNUM - 6 2290 S(IFIELD,IFNUM) = SIGN 2291 GO TO 150 2292C----A NUMBER, FLOAT IT INTO DIGIT 2293 130 DIGIT = IFNUM 2294 ICOL = ICOL + 1 2295C----IS NEXT CHARACTER A SLASH 2296 IF(ICHAR(ICOL).EQ.12) GO TO 135 2297C----NO SLASH, TRANSLATION TERM COMPLETE 2298 ICOL=ICOL-1 2299 GO TO 140 2300C----A SLASH IS NEXT, IS THERE A NUMBER AFTER IT 2301 135 ICOL = ICOL + 1 2302 IFNUM=ICHAR(ICOL) 2303 IF (IFNUM-6) 138,138,9980 2304C----MAKE FRACTION; '5' NOT ALLOWED IN DENOMINATOR 2305 138 IF(IFNUM.EQ.5) GO TO 9980 2306 DIGIT = DIGIT/FLOAT(IFNUM) 2307C----ACCUMULATE THE VECTOR COMPONENT 2308 140 S(IFIELD,4) = S(IFIELD,4) + SIGN*DIGIT 2309C----TEST TO SEE IF THERE ARE MORE CHARACTERS IN THIS SUBFIELD 2310C OR A NEW SUBFIELD 2311 150 ICOL = ICOL + 1 2312 SIGN = 1.0 2313C----TEST IF ALL DONE 2314 IF(ICOL - ICOLMX) 151,151,1000 2315 151 CONTINUE 2316C----A SUBFIELD BEGINS WITH A PLUS OR MINUS 2317 IF(ICHAR(ICOL).EQ.11) GO TO 118 2318 IF(ICHAR(ICOL).EQ.10) GO TO 115 2319C----A COMMA MUST BE NEXT UNLESS ICOL IS ICOLMX+1 WHICH MEANS THE END 2320 IF(ICHAR(ICOL).EQ.13) GO TO 101 2321 163 IF (ICOLMX+1-ICOL) 9980,1000,9980 2322C----EVERYTHING SEEMS OK SEE IF DETERMINANT IS A NICE + OR - 1. 2323 1000 CONTINUE 2324 D=S(1,1)*S(2,2)*S(3,3) 2325 $ -S(1,2)*S(2,1)*S(3,3) 2326 $ -S(1,1)*S(2,3)*S(3,2) 2327 $ -S(1,3)*S(2,2)*S(3,1) 2328 $ +S(1,2)*S(2,3)*S(3,1) 2329 $ +S(1,3)*S(2,1)*S(3,2) 2330 IF(ABS(ABS(D) - 1.0) -.0001) 1001,9985,9985 2331 1001 CONTINUE 2332 IF(IFIELD - 3) 9987,1002,9987 2333 1002 CONTINUE 2334C----LEGAL RETURN 2335 MATSYM = 0 2336 2005 RETURN 2337C 2338C 2339C----ILLEGAL CHARACTER ON SYMMETRY CARD 2340 9975 MATSYM = 1 2341 ICOL = JTXTSC 2342 GO TO 2005 2343C 2344C----SYNTAX ERROR AT NONBLANK CHARACTER ICOL 2345 9980 MATSYM = 2 2346 GO TO 2005 2347C 2348C----DETERMINANT IS NOT + OR - 1. 2349 9985 MATSYM = 4 2350 GO TO 2005 2351C 2352C----INCORRECT NUMBER OF COMMAS 2353 9987 MATSYM = 3 2354 GO TO 2005 2355 END 2356c 2357c 2358c this is subroutine to transfer missetting angle to matrix 2359c ROT = (phiz)*(phiy)*(phix) 2360 subroutine misseting(phi,rot) 2361 real phi(3),rot(3,3) 2362 external cosd, sind 2363 sinx = sind(phi(1)) 2364 cosx = cosd(phi(1)) 2365 siny = sind(phi(2)) 2366 cosy = cosd(phi(2)) 2367 sinz = sind(phi(3)) 2368 cosz = cosd(phi(3)) 2369 rot(1,1) = cosz*cosy 2370 rot(2,1) = sinz*cosy 2371 rot(3,1) = -siny 2372 rot(1,2) = cosz*siny*sinx - sinz*cosx 2373 rot(2,2) = sinz*siny*sinx + cosz*cosx 2374 rot(3,2) = cosy*sinx 2375 rot(1,3) = cosz*siny*cosx + sinz*sinx 2376 rot(2,3) = sinz*siny*cosx - cosz*sinx 2377 rot(3,3) = cosy*cosx 2378 end 2379c if you have a rigid movement x2=R*x1+1 2380c where R(3,3) is the rotation and t(3) is the translation 2381c then x1 =R2*x1-t2 2382c where R2=R-1 2383c t2 =R-1*t 2384c The program reads in R,T matrix (3,4) 1 and 2, then returns an RT 2385c which is rt = rt1*rt2 2386c 2387c 2388c 2389 subroutine morert(rt1,rt2,rt) 2390 real rt(3,4),rt2(3,4),rt1(3,4) 2391 real gt(4,4),gt2(4,4),gt1(4,4) 2392 real extr(4) 2393 real b1(3),b2(3) 2394 integer me(3) 2395c... data statements. Separate declaration and init required for f2c 2396 data extr /0.,0.,0.,1./ 2397c 2398 do i = 1, 4 2399 gt1(i,4) = extr(i) 2400 gt2(i,4) = extr(i) 2401 end do 2402c 2403 do j = 1, 3 2404 do i = 1, 3 2405 gt1(i,j) = rt1(i,j) 2406 gt2(i,j) = rt2(i,j) 2407 end do 2408 end do 2409c 2410 do j = 1, 3 2411 gt1(4,j) = rt1(j,4) 2412 gt2(4,j) = rt2(j,4) 2413 end do 2414c 2415 call matmult(4,4,4,4,gt1,gt2,gt) 2416c 2417 do j = 1, 3 2418 do i = 1, 3 2419 rt(i,j) = rt(i,j) 2420 end do 2421 end do 2422c 2423 do j = 1, 4 2424 rt(j,4) = gt(4,j) 2425 end do 2426c 2427c WRITE(6,*) 'gt' 2428c write(6,'(4f10.5)') gt1 2429c write(6,'(4f10.5)') gt2 2430c write(6,'(4f10.5)') gt 2431c 2432 end 2433c 2434c 2435c a subroutine to transfer a normalized matrix to denzo missetting angle 2436c definition [rot]=[rotx]*[roty]*[rotz] 2437c GuoGuang 930706 2438 subroutine mtodenmis(rot0,phi) 2439 real rot0(3,3) 2440 real rot(3,3),phi(3) 2441 external asind, atand, cosd 2442c sinx = sind(phi(1)) 2443c cosx = cosd(phi(1)) 2444c siny = sind(phi(2)) 2445c cosy = cosd(phi(2)) 2446c sinz = sind(phi(3)) 2447c cosz = cosd(phi(3)) 2448c rot0(1,1) = cosy*cosz 2449c rot0(2,1) = -cosx*sinz+sinx*siny*cosz 2450c rot0(3,1) = sinx*sinz+cosx*siny*cosz 2451c rot0(1,2) = cosy*sinz 2452c rot0(2,2) = cosx*cosz+sinx*siny*sinz 2453c rot0(3,2) = -sinx*cosz+cosx*siny*sinz 2454c rot0(1,3) = -siny 2455c rot0(2,3) = sinx*cosy 2456c rot0(3,3) = cosx*cosy 2457c 2458 call antiarr(3,3,rot0,rot) 2459c sinx = sind(phi(1)) 2460c cosx = cosd(phi(1)) 2461c siny = sind(phi(2)) 2462c cosy = cosd(phi(2)) 2463c sinz = sind(phi(3)) 2464c cosz = cosd(phi(3)) 2465c rot(1,1) = cosz*cosy 2466c rot(2,1) = sinz*cosy 2467c rot(3,1) = -siny 2468c rot(1,2) = cosz*siny*sinx - sinz*cosx 2469c rot(2,2) = sinz*siny*sinx + cosz*cosx 2470c rot(3,2) = cosy*sinx 2471c rot(1,3) = cosz*siny*cosx + sinz*sinx 2472c rot(2,3) = sinz*siny*cosx - cosz*sinx 2473c rot(3,3) = cosy*cosx 2474 phi(2) = asind(-rot(3,1)) 2475 if (abs(phi(2)).ne.90.) then 2476 if (rot(1,1).eq.0.) then 2477 phi(3) = 90. 2478 else 2479 phi(3) = atand(rot(2,1)/rot(1,1)) 2480 if (phi(3).gt.0..and.rot(1,1)/cosd(phi(2)).lt.0.) 2481 1 phi(3) = 180. + phi(3) 2482 if (phi(3).lt.0..and.rot(2,1)/cosd(phi(2)).gt.0.) 2483 1 phi(3) = 180. + phi(3) 2484 end if 2485 if (rot(3,3).eq.0.) then 2486 phi(1) = 90. 2487 else 2488 phi(1) = atand(rot(3,2)/rot(3,3)) 2489 if (phi(1).gt.0..and.rot(3,3)/cosd(phi(2)).lt.0.) 2490 1 phi(1) = 180. + phi(1) 2491 if (phi(1).lt.0..and.rot(3,2)/cosd(phi(2)).gt.0.) 2492 1 phi(1) = 180. + phi(1) 2493 end if 2494 else 2495 WRITE(6,*) 'not implemented yet' 2496 end if 2497 end 2498c 2499c 2500 SUBROUTINE MTOHUBER(E,HB,HB1) 2501 DIMENSION E(3,3),HB(3),HB1(3) 2502C 2503C * ROH ANGLES: -180 < PSI < 180 2504C -90 < TH < 90 2505C -180 < PHI < 180 2506C 2507 ARC=3.14159265359/180.0 2508 if (e(3,2).gt.1.0) e(3,2) = 1.0 2509 if (e(3,2).lt.-1.0) e(3,2) = -1.0 2510210 TH = ASIN(E(3,2)) 2511 COSTH = COS(TH) 2512c IF (COSTH.GE.0.01) GO TO 212 2513 if (costh.lt.1.0e-6) then 2514 ph = 0. 2515 ps = acos(e(1,1)) 2516 if (e(2,1).lt.0) ps = -ps 2517 goto 1500 2518 end if 2519c WRITE (6,1212) COSTH 2520c1212 FORMAT('-*** WARNING : COS(THETA) = ', E8.2, ' ***') 2521212 PH = E(3,3)/COSTH 2522C IF (PH.GT.1.0.OR.PH.LT.-1.0) WRITE(6,1291) PH 2523c IF (ABS(PH).GT.(1.00001)) WRITE(6,1291) PH 25241291 FORMAT(' ARCOS(E33/COSTH)=',F12.5) 25251020 FORMAT(' *** ERROR IN ROH *** ',2F10.4/) 2526 IF (PH.GT.1.0) PH=1.0 2527 IF (PH.LT.-1.0) PH=-1.0 2528 PH=ACOS(PH) 2529 TST=-COSTH*SIN(PH) 2530 IF (TST*E(3,1).LT.0.0) PH=-PH 2531 PS = E(2,2)/COSTH 2532c IF (PS.GT.1.0.OR.PS.LT.-1.0) WRITE(6,1292) PS 2533c if (abs(ps).gt.1.0001) write(6,1292) ps 25341292 FORMAT(' ARCOS(E22/COSTH)=',F12.5) 2535 IF (PS.GT.1.0) PS=1.0 2536 IF (PS.LT.-1.0) PS=-1.0 2537 PS=ACOS(PS) 2538 TST=-COSTH*SIN(PS) 2539 IF (TST*E(1,2).LT.0.0) PS=-PS 2540 TST = -SIN(PS)*SIN(PH)*SIN(TH) + COS(PS)*COS(PH) 2541c IF (ABS(TST-E(1,1)) .GT. 0.1) WRITE (6,1020) TST,E(1,1) 2542 TST = -COS(PS)*SIN(TH)*COS(PH) + SIN(PS)*SIN(PH) 2543c IF (ABS(TST-E(2,3)) .GT. 0.1) WRITE (6,1020) TST,E(2,3) 25441500 continue 2545 PS = PS/ARC 2546 TH = TH/ARC 2547 PH = PH/ARC 2548 PS1 = 180. + PS 2549 TH1 = 180. - TH 2550 PH1 = 180. + PH 2551 CALL KABMOD(PS,TH,PH,PS1,TH1,PH1) 2552 HB(1)=PS 2553 HB(2)=TH 2554 HB(3)=PH 2555 HB1(1)=PS1 2556 HB1(2)=TH1 2557 HB1(3)=PH1 2558 END 2559c 2560c 2561 SUBROUTINE MTOHUBERARC(E,HB,HB1) 2562 DIMENSION E(3,3),HB(3),HB1(3) 2563C 2564C * ROH ANGLES: -180 < PSI < 180 2565C -90 < TH < 90 2566C -180 < PHI < 180 2567C 2568C ARC=3.14159265359/180.0 2569C ARC=1.0 2570210 TH = ASIN(E(3,2)) 2571 COSTH = COS(TH) 2572 IF (COSTH.GE.0.01) GO TO 212 2573 WRITE (6,1212) COSTH 25741212 FORMAT('-*** WARNING : COS(THETA) = ', E8.2, ' ***') 2575212 PH = E(3,3)/COSTH 2576c IF (PH.GT.1.0.OR.PH.LT.-1.0) WRITE(6,1291) PH 25771291 FORMAT(' ARCOS(E33/COSTH)=',F12.5) 25781020 FORMAT(' *** ERROR IN ROH *** ',2F10.4/) 2579 IF (PH.GT.1.0) PH=1.0 2580 IF (PH.LT.-1.0) PH=-1.0 2581 PH=ACOS(PH) 2582 TST=-COSTH*SIN(PH) 2583 IF (TST*E(3,1).LT.0.0) PH=-PH 2584 PS = E(2,2)/COSTH 2585c IF (PS.GT.1.0.OR.PS.LT.-1.0) WRITE(6,1292) PS 25861292 FORMAT(' ARCOS(E22/COSTH)=',F12.5) 2587 IF (PS.GT.1.0) PS=1.0 2588 IF (PS.LT.-1.0) PS=-1.0 2589 PS=ACOS(PS) 2590 TST=-COSTH*SIN(PS) 2591 IF (TST*E(1,2).LT.0.0) PS=-PS 2592 TST = -SIN(PS)*SIN(PH)*SIN(TH) + COS(PS)*COS(PH) 2593c IF (ABS(TST-E(1,1)) .GT. 0.1) WRITE (6,1020) TST,E(1,1) 2594 TST = -COS(PS)*SIN(TH)*COS(PH) + SIN(PS)*SIN(PH) 2595c IF (ABS(TST-E(2,3)) .GT. 0.1) WRITE (6,1020) TST,E(2,3) 2596C PS = PS/ARC 2597C TH = TH/ARC 2598C PH = PH/ARC 2599 PS1 = 180. + PS 2600 TH1 = 180. - TH 2601 PH1 = 180. + PH 2602 CALL KABMOD(PS,TH,PH,PS1,TH1,PH1) 2603 HB(1)=PS 2604 HB(2)=TH 2605 HB(3)=PH 2606 HB1(1)=PS1 2607 HB1(2)=TH1 2608 HB1(3)=PH1 2609 END 2610c 2611c 2612c this is subroutine to transfer missetting angle to matrix 2613c ROT = (phiz)*(phiy)*(phix) 2614 subroutine mtomisset(rot,phi) 2615 real phi(3),rot(3,3) 2616 external asind, atand, cosd 2617c sinx = sind(phi(1)) 2618c cosx = cosd(phi(1)) 2619c siny = sind(phi(2)) 2620c cosy = cosd(phi(2)) 2621c sinz = sind(phi(3)) 2622c cosz = cosd(phi(3)) 2623c rot(1,1) = cosz*cosy 2624c rot(2,1) = sinz*cosy 2625c rot(3,1) = -siny 2626c rot(1,2) = cosz*siny*sinx - sinz*cosx 2627c rot(2,2) = sinz*siny*sinx + cosz*cosx 2628c rot(3,2) = cosy*sinx 2629c rot(1,3) = cosz*siny*cosx + sinz*sinx 2630c rot(2,3) = sinz*siny*cosx - cosz*sinx 2631c rot(3,3) = cosy*cosx 2632 phi(2) = asind(-rot(3,1)) 2633 if (abs(phi(2)).ne.90.) then 2634 if (rot(1,1).eq.0.) then 2635 phi(3) = 90. 2636 else 2637 phi(3) = atand(rot(2,1)/rot(1,1)) 2638 if (phi(3).gt.0..and.rot(1,1)/cosd(phi(2)).lt.0.) 2639 1 phi(3) = 180. + phi(3) 2640 if (phi(3).lt.0..and.rot(2,1)/cosd(phi(2)).gt.0.) 2641 1 phi(3) = 180. + phi(3) 2642 end if 2643 if (rot(3,3).eq.0.) then 2644 phi(1) = 90. 2645 else 2646 phi(1) = atand(rot(3,2)/rot(3,3)) 2647 if (phi(1).gt.0..and.rot(3,3)/cosd(phi(2)).lt.0.) 2648 1 phi(1) = 180. + phi(1) 2649 if (phi(1).lt.0..and.rot(3,2)/cosd(phi(2)).gt.0.) 2650 1 phi(1) = 180. + phi(1) 2651 end if 2652 else 2653 WRITE(6,*) 'not implemented yet' 2654 end if 2655 end 2656c 2657C 2658C * POLAR ANGLES 2659C 2660 SUBROUTINE MTOPOLOR(E,POL,POL1) 2661 DIMENSION E(3,3), POL(3),POL1(3) 2662 ARC=3.14159265359/180.0 2663230 COSPC = (E(1,1)+E(2,2)+E(3,3)-1.0)/2.0 2664 IF (COSPC .LT.-1.0) 2665 1 WRITE(6,1290) 26661290 FORMAT(1X,' SUM OF DIAGONAL-ELEMENTS LESS 2667 1 THAN -1./IS SET TO -1.') 2668 IF (COSPC .LT.-1.0) COSPC=-1.0 2669C --- PC == KAPPA (ALIAS CHI) 2670 PC = ACOS(COSPC) 2671 ARG = (E(2,2)-COSPC)/(1.0-COSPC) 2672 IF (ARG .GE. 0.0) GO TO 235 2673 WRITE(6,1235) ARG 26741235 FORMAT('-*** ARG = ',E10.4, ' * ARG=0 ASSUMED ***') 2675 ARG = 0.0 2676235 COSPA = SQRT(ARG) 2677C --- PA == THETA (ALIAS PSI) 2678 PA = ACOS(COSPA) 2679 SINPA = SIN(PA) 2680 TST = 2*COSPA*SIN(PC) 2681 IF ((E(1,3)-E(3,1))*TST .LT. 0.0) PC=-PC 2682 SINPC = SIN(PC) 2683 IF (SINPA.NE.0.0) GO TO 237 2684 WRITE (6,1237) 26851237 FORMAT('-*** SIN(THETA) = 0; PHI UNDETERMINED ***') 2686 PB = 0. 2687 GO TO 250 2688C --- PB == PHI 2689237 IF(ABS(SINPC) .LT. 0.0001 .AND. 2690 1 ABS(COSPA) .LT. 0.0001) GOTO 245 2691 IF ( ABS(SINPA*SINPC) .LT. 0.001) GO TO 240 2692 PB = ACOS( (E(3,2)-E(2,3)) / (2*SINPA*SINPC) ) 2693 TST = 2.0*SINPA*SIN(PB)*SINPC 2694 IF (TST*(E(1,2)-E(2,1)) .LT. 0.0) PB=-PB 2695 GO TO 250 2696240 FAC = 2.0*SINPA*COSPA*(1.0-COSPC) 2697 PB = ACOS((E(1,2)+E(2,1))/FAC) 2698 TST = FAC*SIN(PB) 2699 IF (-(E(2,3)+E(3,2))*TST .LT. 0.0) PB = -PB 2700 GO TO 250 2701245 DENOM = (1.0-ARG)*(COSPC-1.0) 2702 PB = 0.5*ASIN((E(1,3)+E(3,1))/DENOM) 2703 SINSPB = (COSPC-E(3,3))/DENOM 2704 IF (SIN(PB)**2 .LT. SINSPB) PB = 90.*ARC - PB 2705250 PA = PA/ARC 2706 PB = PB/ARC 2707 PC = PC/ARC 2708 TH1 = 180. - PA 2709 PH1 = 180. + PB 2710 PK1 = - PC 2711 CALL KABMOD(PA,PB,PC,TH1,PH1,PK1) 2712 POL(1)=PA 2713 POL(2)=PB 2714 POL(3)=PC 2715 POL1(1)=TH1 2716 POL1(2)=PH1 2717 POL1(3)=PK1 2718 2719 END 2720c 2721c 2722 subroutine mtopolors(e,p1,p2) 2723 dimension e(3,3),p1(3),p2(3) 2724 dimension b1(3),b2(3) 2725 external acosd 2726 2727 call mtovec(e,b1,p1(3)) 2728C 2729 if (p1(3).eq.0.) then 2730 p1(1) = 0. 2731 p1(2) = 0. 2732 p2(1) = 0. 2733 p2(2) = 0. 2734 p2(3) = 0. 2735 return 2736 end if 2737C 2738 p1(1) = acosd(b1(2)) 2739 dos = sqrt(b1(1)*b1(1)+b1(3)*b1(3)) 2740 if (dos.eq.0.) then 2741 if (vem(3,b1).eq.0.and.p1(3).ne.0.) then 2742 WRITE(6,*) 'vec:',b1,'kapa:',p1(3) 2743 end if 2744 p1(2) = 0. 2745 p2(2) = 0. 2746 p2(1) = 180. - p1(1) 2747 if (p2(1).ge.360.) p2(1) = 0. 2748 return 2749 end if 2750c 2751c 2752 p1(2) = acosd( b1(1)/dos ) 2753c b22 = asind( b1(3)/sqrt(b1(1)*b1(1)+b1(3)*b1(3)) ) 2754 if (b1(3).gt.0.) p1(2) = -p1(2) 2755c 2756 p2(3) = -p1(3) 2757 p2(2) = 180. + p1(2) 2758 if (p2(2).gt.180) p2(2) = p2(2) -360. 2759 p2(1) = 180. - p1(1) 2760 end 2761c 2762c 2763 subroutine mtopolorz(e,p1,p2) 2764c 2765 real e(3,3),p1(3),p2(3) 2766 real b1(3),b2(3) 2767 external acosd 2768 call mtovec(e,b1,p1(3)) 2769c 2770 if (abs(p1(3)).lt.1e-3) then 2771 p1(1) = 0 2772 p2(1) = 0 2773 p1(2) = 0 2774 p2(2) = 0 2775 p2(3) = 0 2776 return 2777 end if 2778 2779c p1(1) = acosd(b1(2)) 2780c p1(2) = acosd( b1(1)/sqrt(b1(1)*b1(1)+b1(3)*b1(3)) ) 2781c b22 = asind( b1(3)/sqrt(b1(1)*b1(1)+b1(3)*b1(3)) ) 2782 p1(1) = acosd(b1(3)) 2783c 2784 if (p1(1).eq.0.or.p1(1).eq.180.) then 2785 p1(2) = 0. 2786 p2(2) = 0. 2787 p2(1) = p1(1) + 180. 2788 if (p2(1).ge.360) p2(1) = p2(1) - 360. 2789 p2(3) = -p1(3) 2790 return 2791 end if 2792c 2793 p1(2) = acosd( b1(1)/sqrt(b1(1)*b1(1)+b1(2)*b1(2)) ) 2794c b22 = asind( b1(2)/sqrt(b1(1)*b1(1)+b1(3)*b1(3)) ) 2795c if (b22.gt.0.) p1(2) = -p1(2) 2796 if (b1(2).lt.0.) p1(2) = -p1(2) 2797c 2798 p2(3) = -p1(3) 2799 p2(1) = 180. - p1(1) 2800 p2(2) = 180. + p1(2) 2801 if (p2(2).gt.180) p2(2) = p2(2) -360. 2802 end 2803c 2804c 2805 SUBROUTINE MTOR_B(IOPT2,E,ROT1,ROT2) 2806C 2807C * R&B ANGLES: -180 < AL < 180 2808C 0 < BE < 180 2809C -180 < GA < 180 2810C 2811 EQUIVALENCE (PS,AL,PA), (TH,BE,PB), (PH,GA,PC) 2812 EQUIVALENCE (S1,SIPS), (S2,SITH), (S3,SIPH) 2813 * ,(C1,COPS), (C2,COTH), (C3,COPH) 2814 DIMENSION E(3,3),E1(3,3),E2(3,3),X(3),X1(3) 2815 DIMENSION ROT1(3),ROT2(3) 2816 DATA E2 / 1., 3*0., 1., 3*0., 1./ 2817 ARC=3.14159265359/180.0 2818 220 BE = ACOS(E(3,3)) 2819 GA = E(3,2)/SIN(BE) 2820c IF (GA.GT.1.0.OR.GA.LT.-1.0) WRITE(6,1293) GA 28211293 FORMAT(' ARCOS(E32/SINBE)=',F12.5) 2822 IF (GA.GT.1.0) GA=1.0 2823 IF (GA.LT.-1.0) GA=-1.0 2824 GA=ACOS(GA) 2825 TST= SIN(BE)*SIN(GA) 2826 IF (TST*E(3,1).LT.0.0) GA=-GA 2827 AL = -E(2,3)/SIN(BE) 2828c IF (AL.GT.1.0.OR.AL.LT.-1.0) WRITE(6,1294) AL 28291294 FORMAT(' ARCOS(-E23/SINBE)=',F12.5) 2830 IF (AL.GT.1.0) AL=1.0 2831 IF (AL.LT.-1.0) AL=-1.0 2832 AL=ACOS(AL) 2833 TST= SIN(BE)*SIN(AL) 2834 IF (TST*E(1,3).LT.0.0) AL=-AL 2835 TST = -SIN(AL)*COS(BE)*SIN(GA) + COS(AL)*COS(GA) 2836 IF (ABS(TST-E(1,1)) .GT. 0.1) WRITE (6,1010) TST,E(1,1) 2837 TST = COS(AL)*COS(BE)*COS(GA) - SIN(AL)*SIN(GA) 2838 IF (ABS(TST-E(2,2)) .GT. 0.1) WRITE (6,1010) TST,E(2,2) 2839 1010 FORMAT(' *** ERROR IN R&B *** ',2F10.4/) 2840 AL = AL/ARC 2841 BE = BE/ARC 2842 GA = GA/ARC 2843 IF (IOPT2.NE.4) GO TO 225 2844 AL=AL-90.0 2845 GA=GA+90.0 2846 225 CONTINUE 2847 AL1 = 180. + AL 2848 BE1 = - BE 2849 GA1 = 180. + GA 2850 CALL KABMOD(AL,BE,GA,AL1,BE1,GA1) 2851C IF (IOPT2.EQ.2) WRITE (6,1220) AL,BE,GA,AL1,BE1,GA1 2852C IF (IOPT2.EQ.4) WRITE (6,1221) AL,BE,GA,AL1,BE1,GA1 2853C GO TO 100 2854 ROT1(1)=AL 2855 ROT1(2)=BE 2856 ROT1(3)=GA 2857 ROT2(1)=AL1 2858 ROT2(2)=BE1 2859 ROT2(3)=GA1 2860 END 2861c 2862c 2863 subroutine mtovec(a,vec,vkapa) 2864c this is a subroutine which give a transfer a orthogonal matrix 2865c to a polar angle rotate a direction of a vector which is 1. 2866c a is the input 3*3 orthogonal matrix. 2867c vec is the output 3-dim vector which indicate the axis direction 2868c kappa is the output polar angle. 2869c Guoguang 901001 bmc 2870 dimension a(3,3),vec(3),a1(3,3),a2(9),vec1(3) 2871 integer*2 sgn(3,8) 2872 external acosd 2873 data sgn/ 1,1,1, 1,1,-1, 1,-1,1, 1,-1,-1, 2874 1 -1,1,1, -1,1,-1, -1,-1,1, -1,-1,-1/ 2875 2876c 2877 trace = a(1,1) + a(2,2) + a(3,3) 2878c write(6,*),' trace ', trace 2879 if (trace .gt.3.) then 2880 if ( (trace-3.) .gt. 1e-4) write(6,*) 'The trace of the '// 2881 1 'matrix is',trace,' more than 3. I set to 3 here.' 2882 trace = 3. 2883 else if (trace .lt.-1.) then 2884 if (( trace+1 ) .lt. -1e-4) write(6,*) 'The trace of the '// 2885 1 'matrix is',trace,' less than -1. I set to -1 here.' 2886 trace = -1. 2887 end if 2888 2889 coskapa = (trace - 1.) / 2. 2890 vkapa = acosd(coskapa) 2891c write(6,*),' coskapa, kapa ', coskapa, kapa 2892c 2893 if (abs(vkapa).lt.0.03) then 2894 vec(1) = 0. 2895 vec(2) = 0. 2896 vec(3) = 0. 2897 return 2898 else 2899 do i = 1, 3 2900 temp = (a(i,i)-coskapa)/(1.-coskapa) 2901 if (temp.lt.0.) then 2902 if (temp.lt.-1e-4.and.(coskapa.ne.-1.)) then 2903 write(6,*) 'Warning: your matrix might be not a orthogonal'// 2904 1 ' one' 2905 write(6,'(a,i1,a,i1,a)') 'vec(',i,')**2 =' 2906 write(6,*) temp,' coskapa,kapa,trace',coskapa,vkapa,trace 2907 end if 2908 temp = 0. 2909 end if 2910 vec1(i) = sqrt(temp) 2911 end do 2912 vv = vem(3,vec1) 2913cc if (abs(vv-1.).gt.1e-4) then 2914c write(6,*) 'Warning: modulus of vector is not 0, but',vv 2915c end if 2916 call arrmc(3,1,vec1,1./vv,vec1) 2917 end if 2918c 2919c if (coskapa.eq.-1.) return 2920 errmin = 1.0 2921 ireal = 0 2922 do is = 1, 8 2923 do i =1, 3 2924 vec(i) = vec1(i) * sgn(i,is) 2925 end do 2926c 2927 call rotvsvec(vec,vkapa,a1) 2928 call arrps(3,3,a1,a,a2) 2929c 2930 err = 0. 2931 do ip = 1,9 2932 err = abs(a2(ip)) + err 2933 end do 2934 if (err.lt.errmin) ireal = is 2935 errmin = min(errmin,err) 2936c 293710 end do 2938c 2939 if (ireal. eq. 0) then 2940c 2941 write(6,*) 'something wrong in this program a,a1' 2942 do is = 1, 8 2943 do i =1, 3 2944 vec(i) = vec1(i) * sgn(i,is) 2945 end do 2946c 2947 call rotvsvec(vec,vkapa,a1) 2948 call arrps(3,3,a1,a,a2) 2949 write(6,*) 'a',a 2950 write(6,*) 'a1',a1 2951 write(6,*) 'a2',a2 2952 err = 0. 2953 do ip = 1,9 2954 err = abs(a2(ip)) + err 2955 end do 2956 write(6,*) vec,vkapa,'vk' 2957 write(6,*) 'Error',err 2958c 2959 end do 2960 stop ' ' 2961 else 2962 do i =1, 3 2963 vec(i) = vec1(i) * sgn(i,ireal) 2964 end do 2965 if (errmin.gt.0.01) then 2966 write(6,*) 'ireal =', ireal 2967 call rotvsvec(vec,vkapa,a1) 2968 call arrps(3,3,a1,a,a2) 2969 write(6,*) a 2970 write(6,*) a1 2971 write(6,*) a2 2972 write(6,*) vec 2973 end if 2974 end if 2975c check 2976c v2 = asind(sinkapa) 2977c if (coskapa.lt.0) v2 = -180. - v2 2978c 2979c WRITE(6,*) 'kapa1,kapa2,error',v2,vkapa,abs(v2-vkapa) 2980c 2981 end 2982c 2983c 2984 subroutine nodir(name1,name2,len2) 2985c a subroutine which extracts the file name from a complete filename. 2986c e.g. input name1 = /nfs/pdb/full/1cnd.pdb 2987c output name2 = 1cnd.pdb 2988c len2 is the length of the second word 2989 character*(*) name1,name2 2990c i1 = index(name1,'.pdb') 2991 i1 = lenstr(name1) 2992 do i = i1, 1, -1 2993 if (name1(i:i).eq.'/') goto 20 2994 end do 2995 i = 0 299620 continue 2997 len2 = i1 - i 2998 name2(1:len2) = name1(i+1:i1) 2999 end 3000c 3001c 3002 SUBROUTINE ORIEN ( NATM, XYZ1, XYZ2, A ) 3003c The subroutine selects three atoms to calculate the initial superposing 3004c matrix A from two molecules in 3*4 array XYZ1 and XYZ2 individually. 3005c Parameter description: 3006c NATM is the number of the atoms 3007c DIMENSION XYZ1(3,NATM), XYZ2(3,NATM), A(3,3) 3008c XYZ1(3,NATM) represent the XYZ of molecule 1 3009c XYZ2(3,NATM) represent the XYZ of molecule 2 3010c A is the output superposing matrix. 3011c COMMON/IAT/ IAT1,IAT2,IAT3,IAT 3012c DATA IAT/0/ 3013c Atom IAT1, IAT2 and IAT3 is used to calculate the matrix. It could be 3014c decided outside the routine if IAT not eq 0. 3015c 3016c 3017 COMMON/IAT/ IAT1,IAT2,IAT3,IAT 3018c DATA IAT/0/,IAT1/1/,IAT2/2/,IAT3/3/ 3019 DIMENSION XYZ1(3,NATM), XYZ2(3,NATM), A(3,3) 3020 DIMENSION V11(3),V12(3) 3021 DIMENSION V21(3),V22(3) 3022 DIMENSION C1(3,3),C2(3,3),C1P(3,3) 3023 DIMENSION B1(3),B2(3) 3024c local elements 3025 REAL T(3) 3026c 3027c If IAT = 0, define the three atoms used to calculate the 3028c initial orientation matrix 3029c 3030 IF (IAT.EQ.0) THEN 3031 IAT1 = 1 3032 IF (NATM.GE.6) THEN 3033 IAT2 = INT(NATM/3) + 1 3034 IAT3 = 2*INT(NATM/3) + 1 3035 ELSE 3036 IAT2 = 2 3037c IAT3 = 3 3038 IAT3 = NATM 3039 END IF 3040 END IF 3041c 304215 if (dist(xyz1(1,iat1),xyz1(1,iat2)).lt.1.e-3.or. 3043 + dist(xyz2(1,iat1),xyz2(1,iat2)).lt.1.e-3) then 3044 iat1=iat1+1 3045 if (iat1.lt.iat2) goto 15 3046 call elize(3,a) 3047 print*, 'Warning: not give any matrix: -orien' 3048 return 3049 end if 305016 if (dist(xyz1(1,iat2),xyz1(1,iat3)).lt.1.e-3.or. 3051 + dist(xyz2(1,iat2),xyz2(1,iat3)).lt.1.e-3) then 3052 iat2=iat2+1 3053 if (iat2.lt.iat3) goto 15 3054 call elize(3,a) 3055 print*, 'Warning: not give any matrix: -orien' 3056 end if 3057C 3058C Calculate the initial matrix and the Eulerian angle. 3059c 3060C 3061C calculate the matrix 1 306220 continue 3063 CALL ARRPS(3,1,XYZ1(1,IAT2),XYZ1(1,IAT1),V11) 3064 CALL ARRPS(3,1,XYZ1(1,IAT2),XYZ1(1,IAT3),V12) 3065 IF (VEM(3,V12).LT.1E-2) GOTO 50 3066 CALL ARRMC (3,1,V12,1./VEM(3,V12),C1(1,1)) 3067 CALL VECCRSMLT (V12,V11,B1) 3068 IF (VEM(3,B1).LT.1E-2) GOTO 50 3069 CALL ARRMC(3,1,B1,1./VEM(3,B1),C1(1,3)) 3070 CALL VECCRSMLT (C1(1,3),C1(1,1),C1(1,2)) 3071 3072C calculate the matrix 2 3073 CALL ARRPS(3,1,XYZ2(1,IAT2),XYZ2(1,IAT1),V21) 3074 CALL ARRPS(3,1,XYZ2(1,IAT2),XYZ2(1,IAT3),V22) 3075 IF (VEM(3,V12).LT.1E-2) GOTO 50 3076 CALL ARRMC (3,1,V22,1./VEM(3,V22),C2(1,1)) 3077 CALL VECCRSMLT (V22,V21,B1) 3078 IF (VEM(3,B1).LT.1E-2) GOTO 50 3079 CALL ARRMC(3,1,B1,1./VEM(3,B1),C2(1,3)) 3080 CALL VECCRSMLT (C2(1,3),C2(1,1),C2(1,2)) 3081 3082c calculate the matrix A = c2 * c1^-1 (x2= A * x1) 3083 call ANTIARR(3,3,C1,C1P) 3084c print*, 'iat1,iat2,iat3',iat1,iat2,iat3 3085 CALL MATMULT(3,3,3,3,C2,C1P,A) 3086C 3087 return 308850 iat2 = iat2 + 1 3089 if (iat2.ne.iat3.and.iat2.le.natm) goto 20 3090 WRITE(6,*) 'warning: cannot give the orientation here' 3091 WRITE(6,*) 'natm,iat,iat1,iat2,iat3',natm,iat,iat1,iat2,iat3 3092 call elize(3,a) 3093 call arrvalue(3,t,0.) 3094 END 3095c 3096c 3097C-----------APPENDIX PROGRAM FOR THE NEWGIN SYSTEM ANGLE 3098C WHEN INOPT1=8 3099 SUBROUTINE OXFORD(ROT,E) 3100 DIMENSION E(3,3),ROT(3) 3101 external cosd, sind 31022000 A1=ROT(1) 3103 A2=ROT(2) 3104 A3=ROT(3) 3105 E(1,1)=COSD(A1)*COSD(A2)*COSD(A3)-SIND(A1) 3106 1*SIND(A3) 3107 E(2,1)=SIND(A1)*COSD(A2)*COSD(A3)+COSD(A1) 3108 1*SIND(A3) 3109 E(3,1)=-(SIND(A2)*COSD(A3)) 3110 E(1,2)=-COSD(A1)*COSD(A2)*SIND(A3)-SIND(A1) 3111 1*COSD(A3) 3112 E(2,2)=-SIND(A1)*COSD(A2)*SIND(A3)+COSD(A1) 3113 1*COSD(A3) 3114 E(3,2)=SIND(A2)*SIND(A3) 3115 E(1,3)=COSD(A1)*SIND(A2) 3116 E(2,3)=SIND(A1)*SIND(A2) 3117 E(3,3)=COSD(A2) 3118 END 3119c 3120c 3121 subroutine packexpnd(cell,natom,xyz,nsym,sym) 3122 parameter(maxatm=25000) 3123 real xyz(3,maxatm),sym(3,4,160) 3124 real cell(6),deor(3,3),orth(3,3) 3125 real cells(6),deors(3,3),orths(3,3) 3126 real sym1(3,4,160) 3127 real tmove0(3),ave(3) 3128 real b1(3),b2(3),b3(3),b4(3),b5(3) 3129 real tmove(3,27) 3130c... data statements. Separate declaration and init required for f2c 3131 data tmove 3132 1 /-1.,-1.,-1., -1.,-1.,0., -1.,-1.,1., 3133 2 -1., 0.,-1., -1., 0.,0., -1., 0.,1., 3134 3 -1., 1.,-1., -1., 1.,0., -1., 1.,1., 3135 4 0.,-1.,-1., -0.,-1.,0., -0.,-1.,1., 3136 5 0., 0.,-1., 0. ,0.,0., 0., 0.,1., 3137 6 0., 1.,-1., 0., 1.,0., 0., 1.,1., 3138 7 1.,-1.,-1., 1.,-1.,0., 1.,-1.,1., 3139 8 1., 0.,-1., 1., 0.,0., 1., 0.,1., 3140 9 1., 1.,-1., 1., 1.,0., 1., 1.,1./ 3141c 3142 nmove = 1 3143 call lgg_crystal(CELL,CELLS,DEOR,ORTH,deors,orths) 3144 call averg(3,natom,xyz,ave) 3145c WRITE(6,*) 'Average xyz',ave 3146 do isym = 1, nsym 3147 call matmult(3,3,3,1,deor,ave,b1) 3148 call matmult(3,3,3,1,sym(1,1,isym),b1,b2) 3149 call arrad(3,1,b2,sym(1,4,isym),b3) 3150 call arrgive(3,b3,b4) 3151c WRITE(6,*) 'b3',b4 3152 do i = 1, 3 3153 call frcinside(b4(i)) 3154 b5(i) = b4(i) - b3(i) 3155 end do 3156c WRITE(6,*) 'b4-b5',b4,b5 3157 call arrad(3,1,b5,sym(1,4,isym),tmove0) 3158c WRITE(6,*) 'tmove0',tmove0 3159 call arrgive(9,sym(1,1,isym),sym1(1,1,nmove)) 3160 call arrgive(3,tmove0,sym1(1,4,nmove)) 3161c write(6,*) 'Operation:',nmove 3162c write(6,11) ((sym1(i,j,nmove),j=1,4),i=1,3) 3163 nmove = nmove + 1 3164 call arrgive(12,sym(1,1,isym),sym1(1,1,nmove)) 3165 do jmove = 1,27 3166 if (jmove.eq.14) goto 40 3167 call arrad(3,1,tmove(1,jmove),tmove0,sym1(1,4,nmove)) 3168 do iatom = 1, natom 3169 call matmult(3,3,3,1,deor,xyz(1,iatom),b1) 3170 call matmult(3,3,3,1,sym1(1,1,nmove),b1,b2) 3171 call arrad(3,1,b2,sym1(1,4,nmove),b3) 3172 if (b3(1).ge.0. .and. b3(1).lt.1. .and. 3173 1 b3(2).ge.0. .and. b3(2).lt.1. .and. 3174 2 b3(3).ge.0. .and. b3(3).lt.1. ) then 3175c write(6,*) 'move-iele-iatom',jmove,iele,iatom 3176c WRITE(6,*) 'xyz-iatom',xyz(1,iatom),xyz(3,iatom),xyz(3,iatom) 3177c WRITE(6,*) 'b1-deor',b1 3178c WRITE(6,*) 'b2-symrot',b2 3179c WRITE(6,*) 'b3-jmove',b3 3180c write(6,*) 'Operation:',nmove 3181c write(6,11) ((sym1(i,j,nmove),j=1,4),i=1,3) 3182 nmove = nmove + 1 3183 call arrgive(12,sym(1,1,isym),sym1(1,1,nmove)) 3184 goto 40 3185 end if 318630 end do 318740 end do 3188 end do 318911 format(4f10.4) 3190c 3191 nsym = nmove - 1 3192c write(6,*) 'Total',nsym,' operations can put the molecule into 3193c 1 the unit cell' 3194 call arrgive(nsym*12,sym1,sym) 3195c 3196 end 3197c 3198c 3199 CHARACTER*16 FUNCTION PLTNAM(IDUM3) 3200c CALL SYSTEM('printenv USER > PLTNAM') 3201c call ccpdpn(55,'PLTNAM','UNKNOWN','F',0,0) 3202c read(55,'(a)') pltnam 3203c do i = 1, 16 3204c if (pltnam(i:i).eq.' ') goto 10 3205c end do 3206cc10 continue 3207c do j = i,16 3208c pltnam(j:j) = ' ' 3209c end do 3210 call getlog(pltnam) 3211 call up(pltnam,16) 3212 END 3213c 3214c 3215C------------STANDARD FUNCTION SUBPROGRAM FOR 3216C----------------DOT PRODUCT OF TWO VECTORS 3217 FUNCTION POIMULT(N1,N2,V1,V2) 3218 DIMENSION V1(N1),V2(N2) 3219 POIMULT=0. 3220 IF (N1.NE.N2) GOTO 90 3221 DO 10 I=1,N1 322210 POIMULT=V1(I)*V2(I)+POIMULT 3223 RETURN 322490 PRINT*,'POIMULT> Two vectors do not have same dimension' 3225 END 3226c 3227c 3228 SUBROUTINE POLOR(POL,E) 3229 DIMENSION E(3,3),POL(3) 3230 EQUIVALENCE (S1,SIPS), (S2,SITH), (S3,SIPH) 3231 1 ,(C1,COPS), (C2,COTH), (C3,COPH) 3232 3233C 3234C * POLAR ANGLE MATRIX (TRANSPOSED ACCORDING TO TOLLIN, ERRATA OF R&B) 3235C 3236 ARC=3.14159265359/180.0 3237 PS=POL(1) 3238 TH=POL(2) 3239 PH=POL(3) 3240 3241 SIPS = SIN(PS*ARC) 3242 COPS = COS(PS*ARC) 3243 SITH = SIN(TH*ARC) 3244 COTH = COS(TH*ARC) 3245 SIPH = SIN(PH*ARC) 3246 COPH = COS(PH*ARC) 3247 3248130 S1SQ = S1*S1 3249 C1SQ = C1*C1 3250 CPC = 1.0 - C3 3251 E(1,1) = C3 + S1SQ*C2*C2*CPC 3252 E(2,1) = S1*C1*C2*CPC - S1*S2*S3 3253 E(3,1) =-S1SQ*C2*S2*CPC - C1*S3 3254 E(1,2) = S1*C1*C2*CPC + S1*S2*S3 3255 E(2,2) = C3 + C1SQ*CPC 3256 E(3,2) =-S1*C1*S2*CPC + S1*C2*S3 3257 E(1,3) =-S1SQ*S2*C2*CPC + C1*S3 3258 E(2,3) =-S1*C1*S2*CPC - S1*C2*S3 3259 E(3,3) = C3 + S1SQ*S2*S2*CPC 3260 END 3261c 3262c 3263 subroutine polors(pol,e) 3264c my program to calculate matrix from polar angle 3265c pol is psi phi kappa 3266c where psi is angle between axis and Y 3267c phi is angle between image of axis in XZ and X 3268c kappa is rotating angle 3269 3270 dimension e(3,3),pol(3) 3271 dimension b1(3) 3272 external cosd, sind 3273 b1(2) = cosd(pol(1)) 3274 b1(1) = sind(pol(1))*cosd(pol(2)) 3275 b1(3) = - sind(pol(1))*sind(pol(2)) 3276 call rotvsvec(b1,pol(3),e) 3277 end 3278c 3279c 3280 subroutine polorz(pol,e) 3281c my program to calculate matrix from polar angle (z system) 3282c pol is psi phi kappa 3283c where psi is angle between axis and Z 3284c phi is angle between image of axis in XY and X 3285c kappa is rotating angle 3286 3287 real e(3,3),pol(3) 3288 real b1(3) 3289 external cosd, sind 3290c b1(2) = cosd(pol(1)) 3291c b1(1) = sind(pol(1))*cosd(pol(2)) 3292c b1(3) = - sind(pol(1))*sind(pol(2)) 3293 b1(3) = cosd(pol(1)) 3294 b1(1) = sind(pol(1))*cosd(pol(2)) 3295 b1(2) = sind(pol(1))*sind(pol(2)) 3296 call rotvsvec(b1,pol(3),e) 3297 end 3298c 3299c 3300 SUBROUTINE POS2VEC(NATM,POS,VEC) 3301C This subroutine is to calculate the center-atom vector from the 3302c coordinate. NATM vector will be produced from NATM atoms when NATM >= 3. 3303c The NATM vectors is CEN->atom2, CEN->atom3,....CEN->atomN, where 3304C CEN is the centre of the NATM atoms It is used to superimpose the 3305c orientation of two molecules. 3306c written by Guoguang Lu 3307c 3308C DIMENSION POS(3,NATM), VEC(3,NATM) 3309C Parameter description: 3310C NATM is the number of atoms 3311c POS(1->3,i) is the coordinate of the i'th atom 3312c VEC is the output vectors 3313c 3314 DIMENSION POS(3,NATM), VEC(3,NATM),C(3) 3315 IF (NATM.LT.3) STOP 'Number of the atoms is less than 3' 3316 CALL AVERG(3,NATM,POS,C) 3317C C is the center coordinate of the atoms 3318 DO I=1,NATM 3319 CALL ARRPS(3,1,POS(1,I),C,VEC(1,I)) 3320C CALL ARRMC (3,1,VEC(1,I),1./VEM(3,VEC(1,I)),VEC(1,I)) 3321 END DO 3322 END 3323c 3324C 3325C * R&B - MATRIX 3326C 3327 SUBROUTINE R_B(ROT,E,IOTP) 3328 DIMENSION E(3,3),ROT(3) 3329 EQUIVALENCE (PS,AL,PA), (TH,BE,PB), (PH,GA,PC) 3330 EQUIVALENCE (S1,SIPS), (S2,SITH), (S3,SIPH) 3331 1 ,(C1,COPS), (C2,COTH), (C3,COPH) 3332 3333 PS=ROT(1) 3334 TH=ROT(2) 3335 PH=ROT(3) 3336 3337 IF (IOTP.EQ.4) THEN 3338 PS=PS+90.0 3339 PH=PH-90.0 3340 END IF 3341 3342 ARC=3.14159265359/180.0 3343 SIPS = SIN(PS*ARC) 3344 COPS = COS(PS*ARC) 3345 SITH = SIN(TH*ARC) 3346 COTH = COS(TH*ARC) 3347 SIPH = SIN(PH*ARC) 3348 COPH = COS(PH*ARC) 3349 3350120 E(1,1) =-S1*C2*S3 + C1*C3 3351 E(2,1) = C1*C2*S3 + S1*C3 3352 E(3,1) = S2*S3 3353 E(1,2) =-S1*C2*C3 - C1*S3 3354 E(2,2) = C1*C2*C3 - S1*S3 3355 E(3,2) = S2*C3 3356 E(1,3) = S1*S2 3357 E(2,3) = - C1*S2 3358 E(3,3) = C2 3359 END 3360c 3361c 3362 real function lgg_radii(natm,xyz) 3363c pjx: renamed from radii to avoid clash with common block 3364c of same name in refmac4 3365c Explicitly typed to real 3366c 3367c--calculate average radii of a group of atoms. 3368c first calculate the center of gravity 3369c then the mean distance from each atom to this center. 3370c 3371 real xyz(3,natm) 3372 real b1(3),b2(3) 3373c 3374 lgg_radii = 0. 3375 call averg(3,natm,xyz,b1) 3376 do i = 1, natm 3377 call arrps(3,1,xyz(1,i),b1,b2) 3378 lgg_radii = lgg_radii + vem(3,b2) 3379 end do 3380 lgg_radii = lgg_radii/float(natm) 3381 end 3382c 3383c 3384c this is a subroutine to change orthogonalization matrix according to 3385c cell dimension raxis spindle and xray axis 3386c input 3387c character cell(6) ------ cell dimensions 3388c character aspin*3 ------ raxis spindle axis 3389c character axray*2 ------ raxis xray axis 3390c output 3391c CELLS*6 --- cell dimensions in reciprocal space 3392c DEOR3*3 --- Deorthogonalization matrix in real space 3393c ORTH3*3 --- Orthogonalization matrix in receprocal space. 3394c DEORS3*3 --- Deorthogonalization matrix in reciprocal space 3395c ORTHS3*3 --- Orthogonalization matrix in reciprocal space. 3396c AND COMMON/DET/ DET 3*3 MATRIX transfer from statnd X-a orthogonal matrix 3397c 3398c Guoguang 930723 3399c 3400c PJX 010125 renamed common block det to lggdet to avoid 3401c clash with subroutine det in rsps. 3402c This common block doesn't appear to be used anywhere else? 3403c 3404 subroutine raxcrystl(cell,aspin,axray,cells,deor,orth,deors,orths) 3405 common/lggdet/ det 3406 real det(3,3),buf(3,3),buf2(3,3) 3407 real cell(6),cells(6) 3408 real deor(3,3),orth(3,3) 3409 real deors(3,3),orths(3,3) 3410 real pdir(3,6) 3411 INTEGER IC(6) 3412 real cof(6) 3413 character*3 aspin,pspin(6) 3414 character*2 axray,pxray(6) 3415 DIMENSION B1(3),B2(3),IUF(3) 3416c... data statements. Separate declaration and init required for f2c 3417 data pdir /1.,0.,0.,0.,1.,0.,0.,0.,1., 3418 1 -1.,0.,0.,0.,-1.,0.,0.,0.,-1./ 3419 data ic /1,2,3,1,2,3/ 3420 data cof /1.,1.,1.,-1.,-1.,-1./ 3421 data pspin /'+a*','+b*','+c*','-a*','-b*','-c*'/ 3422 data pxray /'+a','+b','+c','-a','-b','-c'/ 3423c 3424c standard orthogonalization matrix 3425c 3426 call lgg_crystal(cell,cells,deor,orth,deors,orths) 3427c 3428c transferring matrix from new system to old det 3429c 3430c xray axis 3431 do i = 1, 6 3432 if (axray.eq.pxray(i)) then 3433 call arrmc(3,1,orth(1,ic(i)), 3434 1 cof(i)/vem(3,orth(1,ic(i))),det(1,1)) 3435 goto 30 3436 end if 3437 end do 3438 WRITE(6,*) 'X-ray axis',(pxray(i),' ',i=1,6) 3439 WRITE(6,*) 'but you have ',axray 3440 stop 'wrong xray axis' 344130 continue 3442c spindle axis 3443 do i = 1, 6 3444 if (aspin.eq.pspin(i)) then 3445 call arrmc(3,1,orths(1,ic(i)), 3446 1 cof(i)/vem(3,orths(1,ic(i))),det(1,3)) 3447 goto 40 3448 end if 3449 end do 3450 WRITE(6,*) 'spindle axis',(pspin(i),' ',i=1,6) 3451 WRITE(6,*) 'but you have ',aspin 3452 stop 'wrong xray axis' 3453c 345440 continue 3455 call veccrsmlt(det(1,3),det(1,1),det(1,2)) 3456c 3457 call arrgive(9,orth,buf) 3458 call antiarr(3,3,det,buf2) 3459 call matmult(3,3,3,3,buf2,buf,orth) 3460c 3461 CALL ARRGIVE(9,ORTH,DEOR) 3462 CALL IVSN(3,DEOR,B1,B2,IUF,DE,1.0E-6) 3463 CALL ANTIARR(3,3,ORTH,DEORS) 3464 CALL ANTIARR(3,3,DEOR,ORTHS) 3465C 3466c DO I = 1, 3 3467c CELLS(I) = VEM(3,orths(1,I)) 3468c END DO 3469C 3470c CELLS(4) = ANGLE(ORTHS(1,2),ORTHS(1,3)) 3471c CELLS(5) = ANGLE(ORTHS(1,3),ORTHS(1,1)) 3472c CELLS(6) = ANGLE(ORTHS(1,1),ORTHS(1,2)) 3473c 3474 end 3475c 3476c 3477c this is a subroutine to change raxis spindle and xray axis into U matrix 3478c input 3479c character aspin*3 ------ raxis spindle axis 3480c character axray*2 ------ raxis xray axis 3481c real*4 umat(3,3) ------ raxis U matrix 3482c Guoguang 930723 3483c 3484 subroutine raxumat(aspin,axray,umat) 3485 real umat(3,3) 3486 real rmat(3,3) 3487 real pdir(3,6) 3488 character*3 aspin,pspin(6) 3489 character*2 axray,pxray(6) 3490c... data statements. Separate declaration and init required for f2c 3491 data pdir /1.,0.,0.,0.,1.,0.,0.,0.,1., 3492 1 -1.,0.,0.,0.,-1.,0.,0.,0.,-1./ 3493 data pspin /'+a*','+b*','+c*','-a*','-b*','-c*'/ 3494 data pxray /'+a','+b','+c','-a','-b','-c'/ 3495c 3496 do i = 1, 6 3497 if (aspin.eq.pspin(i)) then 3498 call arrgive(3,pdir(1,i),umat(1,3)) 3499 goto 20 3500 end if 3501 end do 3502 WRITE(6,*) ' allowed Spindle axis',(pspin(i),' ',i=1,6) 3503 stop 'wrong spindle axis' 350420 continue 3505c 3506 do i = 1, 6 3507 if (axray.eq.pxray(i)) then 3508 call arrgive(3,pdir(1,i),umat(1,1)) 3509 goto 30 3510 end if 3511 end do 3512 WRITE(6,*) 'X-ray axis',(pxray(i),' ',i=1,6) 3513 stop 'wrong xray axis' 351430 continue 3515 call veccrsmlt(umat(1,3),umat(1,1),umat(1,2)) 3516 call antiarr(3,3,umat,rmat) 3517 call arrgive(9,rmat,umat) 3518c 3519 end 3520c 3521c 3522 SUBROUTINE RECEPICAL(CELL,DEOR,ORTH) 3523C This is a routine which inputs cell parameters, CELL and converts 3524c them into reciprocal cell parameters. It then gives a deorthogonalization 3525c and orthogonalization matrix in reciprocal space. 3526c 3527 DIMENSION DEOR(3,3),ORTH(3,3),CELL(6) 3528 DIMENSION BUFF(3,3) 3529 DIMENSION B1(3),B2(3),IUF(3) 3530 external acosd, cosd, sind 3531 AP=CELL(4) 3532 BA=CELL(5) 3533 GA=CELL(6) 3534 CALL ARRMC(3,3,BUFF,0.,BUFF) 3535 COASTAR=(COSD(BA)*COSD(GA)-COSD(AP))/(SIND(BA) 3536 1*SIND(GA)) 3537 SIASTAR=SQRT(1.-COASTAR*COASTAR) 3538 BUFF(1,1)=1. 3539 BUFF(1,2)=COSD(GA) 3540 BUFF(2,2)=SIND(GA) 3541 BUFF(1,3)=COSD(BA) 3542 BUFF(2,3)=-SIND(BA)*COASTAR 3543 BUFF(3,3)=SIND(BA)*SIASTAR 3544C 3545 DO I =1, 3 3546 CALL ARRMC(3,1,BUFF(1,I),CELL(I),BUFF(1,I)) 3547 END DO 3548c 3549c WRITE(6,*) buff 3550C 3551 VOLUM = VLDIM3(BUFF) 3552C 3553 CALL VECCRSMLT(BUFF(1,2),BUFF(1,3),DEOR(1,1)) 3554 CALL VECCRSMLT(BUFF(1,3),BUFF(1,1),DEOR(1,2)) 3555 CALL VECCRSMLT(BUFF(1,1),BUFF(1,2),DEOR(1,3)) 3556 CALL ARRMC(3,3,DEOR,1./VOLUM,DEOR) 3557C 3558c WRITE(6,*) volum 3559c WRITE(6,*) i,cell(i),vem(3,deor(1,i)) 3560 DO I = 1, 3 3561 CELL(I) = VEM(3,DEOR(1,I)) 3562c WRITE(6,*) CELL(I) 3563c WRITE(6,*) DEOR 3564 CALL ARRMC(3,1,DEOR(1,I),1./CELL(I),DEOR(1,i)) 3565 END DO 3566C 3567 CELL(4) = ACOSD(POIMULT(3,3,DEOR(1,2),DEOR(1,3))) 3568 CELL(5) = ACOSD(POIMULT(3,3,DEOR(1,3),DEOR(1,1))) 3569 CELL(6) = ACOSD(POIMULT(3,3,DEOR(1,1),DEOR(1,2))) 3570 3571C 3572 CALL ARRMC(3,3,DEOR,1.,ORTH) 3573 CALL IVSN(3,ORTH,B1,B2,IUF,VAL,1E-6) 3574 RETURN 3575 END 3576c 3577c 3578 SUBROUTINE RECEPICAL0(CELL,DEOR,ORTH) 3579C This is a routine which inputs cell parameters, CELL and converts 3580c them into reciprocal cell parameters. It then gives a deorthogonalization 3581c and orthogonalization matrix in reciprocal space. 3582c 3583 DIMENSION DEOR(3,3),ORTH(3,3),CELL(6) 3584 DIMENSION BUFF(3,3) 3585 DIMENSION B1(3),B2(3),IUF(3) 3586 external acosd, cosd, sind 3587 AP=CELL(4) 3588 BA=CELL(5) 3589 GA=CELL(6) 3590 CALL ARRVALUE(3*3,BUFF,0.) 3591 COASTAR=(COSD(BA)*COSD(GA)-COSD(AP))/(SIND(BA) 3592 1*SIND(GA)) 3593 SIASTAR=SQRT(1.-COASTAR*COASTAR) 3594 BUFF(1,1)=1. 3595 BUFF(1,2)=COSD(GA) 3596 BUFF(2,2)=SIND(GA) 3597 BUFF(1,3)=COSD(BA) 3598 BUFF(2,3)=-SIND(BA)*COASTAR 3599 BUFF(3,3)=SIND(BA)*SIASTAR 3600C 3601 DO I =1, 3 3602 CALL ARRMC(3,1,BUFF(1,I),CELL(I),BUFF(1,I)) 3603 END DO 3604c 3605c WRITE(6,*) buff 3606C 3607 VOLUM = VLDIM3(BUFF) 3608C 3609 CALL VECCRSMLT(BUFF(1,2),BUFF(1,3),ORTH(1,1)) 3610 CALL VECCRSMLT(BUFF(1,3),BUFF(1,1),ORTH(1,2)) 3611 CALL VECCRSMLT(BUFF(1,1),BUFF(1,2),ORTH(1,3)) 3612 CALL ARRMC(3,3,ORTH,1./VOLUM,ORTH) 3613C 3614c WRITE(6,*) volum 3615c WRITE(6,*) i,cell(i),vem(3,ORTH(1,i)) 3616 DO I = 1, 3 3617 CELL(I) = VEM(3,ORTH(1,I)) 3618c WRITE(6,*) CELL(I) 3619c WRITE(6,*) ORTH 3620 CALL ARRMC(3,1,ORTH(1,I),1./CELL(I),DEOR(1,i)) 3621 END DO 3622C 3623 CELL(4) = ACOSD(POIMULT(3,3,DEOR(1,2),DEOR(1,3))) 3624 CELL(5) = ACOSD(POIMULT(3,3,DEOR(1,3),DEOR(1,1))) 3625 CELL(6) = ACOSD(POIMULT(3,3,DEOR(1,1),DEOR(1,2))) 3626 3627C 3628 CALL ARRGIVE(3*3,ORTH,DEOR) 3629 CALL IVSN(3,DEOR,B1,B2,IUF,VAL,1E-6) 3630 RETURN 3631 END 3632c 3633c 3634 SUBROUTINE REDSTRIN(NUNIT,NCHA,TXT,NPAR) 3635C A subroutine to write a string of characters in unit NUNIT according 3636c to the spaces in the string. It is used to change a character into 3637c parameters under help of a file NUNIT. 3638c NUNIT is the unit number. 3639c NCHA is length of character string. 3640c TXT is the NCHA bytes character. 3641C NPAR is how many parameters in this string. 3642 CHARACTER*120 TXT 3643 NPAR = 0 3644 REWIND (NUNIT) 3645 I = 1 364610 IF (I.LE.NCHA) THEN 3647 IF (TXT(I:I).NE.' ') THEN 3648 J = I 364920 IF (J.GE.NCHA) THEN 3650 WRITE(NUNIT,'(A)') TXT(I:J) 3651 NPAR = NPAR + 1 3652 ELSE IF (TXT(J+1:J+1).EQ.' ') THEN 3653 WRITE(NUNIT,'(A)') TXT(I:J) 3654 NPAR = NPAR + 1 3655 ELSE 3656 J = J + 1 3657 GOTO 20 3658 END IF 3659 I = J + 1 3660 GOTO 10 3661 ELSE IF (I.LT.NCHA) THEN 3662 I = I + 1 3663 GOTO 10 3664 END IF 3665 END IF 3666C 3667 REWIND(NUNIT) 3668 END 3669c 3670c 3671 subroutine refmtoroh(amat,roh,roh2,rms) 3672c This is a subroutine to calculate from a matrix to Eulerian angles in 3673c Rober O Huber system and refine ROH angles to minimize 3674c Sigma((Aroh(i,j)-Amat(i,j)**2). The output will include two 3675c symmetric ROH angles a rms value of input matrix and the matrix 3676c from ROH angles. 3677c Input: amat(3,3) --- Input matrix 3678c Output: roh(3,3) --- output ROH angle 3679c roh2(3,3) --- output ROH angle 3680c rms --- rms value between amat and mat from ROH angle. 3681c Guoguang Lu, 950831. Purdue Univ 3682c 3683 real amat(3,3),roh(3),roh2(3) 3684 real bmat(3,3),b1(3),b2(3),cmat(9) 3685 real rohnew(3) 3686 real vt(9,3),vvv(3,3),delta(3),det(3) 3687c real droh(3,3,3) 3688c 3689 ncyc = 0 3690 deg = 180./3.141592654 3691 call mtohuber(amat,roh,roh2) 3692 call huber(roh,bmat) 3693 call arrps(9,1,bmat,amat,cmat) 3694c WRITE(6,*) 'cmat',cmat 3695 rms = vem(9,cmat) 3696c11 format(3f10.7) 3697c write(6,*) 'bmat' 3698c write(6,11) bmat 3699c 370010 continue 3701c WRITE(6,*) 'cycle,rms',rms,ncyc 3702c WRITE(6,*) 'roh',roh 3703 ncyc = ncyc + 1 3704c 3705 call drvrohthd(1,roh,vt(1,1)) 3706 call drvrohthd(2,roh,vt(1,2)) 3707 call drvrohthd(3,roh,vt(1,3)) 3708c 3709c write(6,*) 'vt' 3710c write(6,'(3f12.7,2x,f12.7)') ((vt(i,j),j=1,3),cmat(i),i=1,9) 3711 call lsqeq(9,3,vt,cmat,det,vvv,b1) 3712c WRITE(6,*) 'det',det 3713 call arrps(3,1,roh,det,rohnew) 3714c WRITE(6,*) 'rohnew',rohnew 3715 call huber(rohnew,bmat) 3716c write(6,*) 'bmatnew' 3717c write(6,11) bmat 3718 call arrps(9,1,bmat,amat,cmat) 3719c WRITE(6,*) 'cmatnew',cmat 3720 rms1 = vem(9,cmat) 3721 if (rms1.lt.rms) then 3722 call arrgive(3,rohnew,roh) 3723 rms = rms1 3724 goto 10 3725 else 3726c WRITE(6,*) 'Refinement finished with rms,rms1:',rms,rms1 3727 roh2(1) = 180. + roh(1) 3728 roh2(2) = 180. - roh(2) 3729 roh2(3) = 180. + roh(3) 3730 end if 3731c 3732 end 3733c 3734c 3735 SUBROUTINE REFORN ( NATM, XYZ1, XYZ2, A, VT, DIS ) 3736C DIMENSION XYZ1(3,NATM), XYZ2(3,NATM), A(3,3) 3737C DIMENSION VT(3*NATM,3), DIS(NATM,3) 3738C A subroutine to give the superimpose matrix and vector from MOL1 to 3739c MOL2 i.e. XYZ2(3*NATM) = A(3*3) * XYZ1(3*NATM) + T(3) 3740C Where NATM represents the number of the atoms in each molecule 3741c which will be superimposed. 3742c XYZ1(3,NATM) represents the coordinates of MOL1 3743c XYZ2(3,NATM) represents the coordinates of MOL2 3744C A(3*3) represents the input initial matrix and 3745c output matrix. 3746C VT(3*6*NATM) is a 6*NATM-dimensional buffer array 3747C DIS(3*NATM) is a 3*NATM-dimensional buffer array 3748c 3749C COMMON/SHIFS/ SCALR,SCALT 3750C Note: There should be a scale of the shifts. i.e. Shift= scale * shifts 3751c Here SCALR is the rotation shiftscale 3752c Here SCALT is the translation shiftscale 3753c The initial and suggested SCALR is 1. 3754c The initial and suggested SCALt is 1. 3755C 3756C COMMON/RMS/ RMS,SMEAN,NREF,NREF1 3757C RMS is the final R.M.S between two molecules. 3758C SMEAN is the mean distance. 3759c NREF is the number of cycles of rotation refinement 3760c NREF1 is the number of cycles of translation refinement 3761c They could be used to judge the SCALR and SCALS 3762c 3763c The routine uses the atom-atom vector to perform the Eulerian angle least 3764c square refinement. Testing shows that that this method might be more 3765c accurate than others at least in orientation. 3766c 3767C SUPOS1 is almost same with SUPOS, but used as a subroutine as a first 3768c step of the refinement by subroutine SUPRIMP. 3769c 3770c by Guoguang Lu 3771cv 3772c added an upper limit for number of refinement cycles to do. 3773c C. Vonrhein 3774cv 3775C 30/01/1989 3776C 3777C 3778 COMMON/RMS/ RMS,SMEAN,NREF,NREF1 3779c COMMON/SHIFS/ SCALR,SCALS 3780 PARAMETER (NATM0 = 50000) 3781c NATM0 is the largest number of atoms. 3782 DIMENSION XYZ1(3,NATM), XYZ2(3,NATM), A(3,3) 3783 DIMENSION VT(3*NATM,3), DIS(3,NATM) 3784 DIMENSION VEC1(3,NATM0),VEC2(3,NATM0),VEC1P(3,NATM0) 3785C 3786 DIMENSION DA(3,3,3),VVV(3,3) 3787 DIMENSION B1(3),B2(3),ROHOLD(3),ROH(3),DETAROH(3),TOLD(3) 3788C MAXREF = maximum number of refinements to do 3789C 3790 INTEGER MAXREF 3791 PARAMETER (MAXREF=100) 3792C 3793c... data statements. Separate declaration and init required for f2c 3794 DATA SCALR/1./,SCALT/1./ 3795C 3796 DEG=180./3.14159265359 3797 NREF=0 3798c NREF1=0 3799C 3800 CALL MTOHUBER(A,ROH,B1) 3801c 3802c calculate the atom-atom vector 3803c 3804 CALL POS2VEC(NATM,XYZ1,VEC1) 3805 CALL POS2VEC(NATM,XYZ2,VEC2) 3806C 3807 CALL MATMULT(3,3,3,NATM,A,VEC1,VEC1P) 3808 CALL ARRPS(3,NATM,VEC1P,VEC2,DIS) 3809C A*V1 - V2 3810 RMS1= DOSQ( 3*NATM, DIS ) 3811 RMS = sqrt(RMS1)/float(NATM) 3812C 381350 CONTINUE 3814C 3815C Refine the superimposing Eulerian Angle' 3816c 3817C WRITE(6,'(A)') '0' 3818C WRITE(6,*) 'Cycle of Eulerian angle refinement ',NREF 3819C WRITE(6,*) 'Superimposing matrix' 3820C WRITE(6,10) ((A(I,J),J=1,3),I=1,3) 3821C WRITE(6,*) 'Eulerian Angle', ROH 3822C10 format (3F10.5) 3823C WRITE(6,*) 'SIGMA (Deta vector)^2 =', RMS 3824C 3825C 3826c 3827C COMPUTE THE DETATH1,DETATH2,DETATH3 than the normal matrix 3828C 3829 DO I = 1, 3 3830 CALL DRVROHTH(I,ROH,DA(1,1,I)) 3831 CALL MATMULT(3,3,3,NATM,DA(1,1,I),VEC1,VT(1,I)) 3832C (DERIV(A)/DERIV(THETAi))*VEC1 3833 END DO 3834c 3835 CALL LSQEQ(NATM*3,3,VT,DIS,DETAROH,VVV,B1) 3836C To change it into degree from arc unit 3837 CALL ARRMC(3,1,DETAROH,DEG,DETAROH) 3838C 3839 CALL ARRMC(3,1,DETAROH,-SCALR,DETAROH) 3840C WRITE(6,*) 'DetaROH = ',DETAROH 3841 CALL ARRMC(3,1,ROH,1.,ROHOLD) 3842C SAVE THE PREVIOUS ROH 3843 CALL ARRAD(3,1,ROHOLD,DETAROH,ROH) 3844C ROH=ROHOLD+DETAROH 3845 CALL HUBER(ROH,A) 3846C NEW MATRIX 3847 CALL MATMULT(3,3,3,NATM,A,VEC1,VEC1P) 3848 CALL ARRPS(3,NATM,VEC1P,VEC2,DIS) 3849C A*V1 - V2 3850 RMS=dosq(3*NATM,DIS) 3851C 3852C do only a maximum of MAXREF refinement cycles 3853C 3854 IF ((NREF.LE.MAXREF).AND.(RMS.LT.RMS1)) THEN 3855c 3856 RMS1=RMS 3857 NREF=NREF+1 3858 GOTO 50 3859 END IF 3860 3861C WRITE(6,'(A)') '0' 3862C WRITE(6,*) 'That is the final rotation result.' 3863C WRITE(6,'(A)') '0' 3864 CALL ARRMC(3,1,ROHOLD,1.,ROH) 3865 CALL HUBER(ROH,A) 3866 3867 END 3868c 3869c 3870 SUBROUTINE REFORNFIN ( NATM, XYZ1, XYZ2, A, VT, DIS ) 3871C DIMENSION XYZ1(3,NATM), XYZ2(3,NATM), A(3,3) 3872C DIMENSION VT(3*NATM,3), DIS(NATM,3) 3873C A subroutine to give the superimpose matrix and vector from MOL1 to 3874c MOL2 i.e. XYZ2(3*NATM) = A(3*3) * XYZ1(3*NATM) + T(3) 3875C Where NATM represents the number of the atoms in each molecule 3876c which will be superimposed. 3877c XYZ1(3,NATM) represents the coordinates of MOL1 3878c XYZ2(3,NATM) represents the coordinates of MOL2 3879C A(3*3) represents the input initial matrix and 3880c output matrix. 3881C VT(3*6*NATM) is a 6*NATM-dimensional buffer array 3882C DIS(3*NATM) is a 3*NATM-dimensional buffer array 3883c 3884C COMMON/SHIFS/ SCALR,SCALT 3885c DATA SCALR/1./,SCALT/1./ 3886C Note: There should be a scale of the shifts. i.e. Shift= scale * shifts 3887c Here SCALR is the rotation shiftscale 3888c Here SCALT is the translation shiftscale 3889c The initial and suggested SCALR is 1. 3890c The initial and suggested SCALt is 1. 3891C 3892C COMMON/RMS/ RMS,SMEAN,NREF,NREF1 3893C RMS is the final R.M.S between two molecules. 3894C SMEAN is the mean distance. 3895c NREF is the number of cycles of rotation refinement 3896c NREF1 is the number of cycles of translation refinement 3897c They could be used to judge the SCALR and SCALS 3898c 3899c The program use the atom-atom vector to perform the Eulerian angle least 3900c square refinement. Testing shows that this method might be more 3901c accurate than others at least in orientation. 3902c 3903C SUPOS1 is almost the same as SUPOS, but used as a subroutine as a first 3904c step of the refinement by subroutine SUPRIMP. 3905c 3906c by Guoguang Lu 3907C 30/01/1989 3908C 3909C 3910 COMMON/RMS/ RMS,SMEAN,NREF,NREF1 3911C COMMON/SHIFS/ SCALR,SCALS 3912 PARAMETER (NATM0 = 50000) 3913c NATM0 is the largest number of atoms. 3914 DIMENSION XYZ1(3,NATM), XYZ2(3,NATM), A(3,3) 3915 DIMENSION VT(3*NATM,3), DIS(3,NATM) 3916 DIMENSION VEC1(3,NATM0),VEC2(3,NATM0),VEC1P(3,NATM0) 3917C 3918 DIMENSION DA(3,3,3),VVV(3,3) 3919 DIMENSION B1(3),B2(3),ROHOLD(3),ROH(3),DETAROH(3),TOLD(3) 3920C 3921C 3922C DEG=180./3.14159265359 3923C 3924 SCALR = 1. 3925 NREF=0 3926c NREF1=0 3927C 3928 CALL MTOHUBERARC(A,ROH,B1) 3929c 3930c calculate the atom-atom vector 3931c 3932 CALL POS2VEC(NATM,XYZ1,VEC1) 3933 CALL POS2VEC(NATM,XYZ2,VEC2) 3934C 3935 CALL MATMULT(3,3,3,NATM,A,VEC1,VEC1P) 3936 CALL ARRPS(3,NATM,VEC1P,VEC2,DIS) 3937C A*V1 - V2 3938 RMS1= DOSQ( 3*NATM, DIS ) 3939 RMS = RMS1 3940C 394150 CONTINUE 3942C 3943C Refine the superimposing Eulerian Angle' 3944c 3945C WRITE(6,'(A)') '0' 3946C WRITE(6,*) 'Cycle of Eulerian angle refinement ',NREF 3947C WRITE(6,*) 'Superimposing matrix' 3948C WRITE(6,10) ((A(I,J),J=1,3),I=1,3) 3949C WRITE(6,*) 'Eulerian Angle', ROH 3950C10 format (3F10.5) 3951C WRITE(6,*) 'SIGMA (Deta vector)^2 =', RMS 3952C 3953C 3954c 3955C COMPUTE THE DETATH1,DETATH2,DETATH3 than the normal matrix 3956C 3957 DO I = 1, 3 3958 CALL DRVROHTHARC(I,ROH,DA(1,1,I)) 3959 CALL MATMULT(3,3,3,NATM,DA(1,1,I),VEC1,VT(1,I)) 3960C (DERIV(A)/DERIV(THETAi))*VEC1 3961 END DO 3962c 3963 CALL LSQEQ(NATM*3,3,VT,DIS,DETAROH,VVV,B1) 3964C To change it into degree from arc unit 3965C CALL ARRMC(3,1,DETAROH,DEG,DETAROH) 3966C 396730 CONTINUE 3968 CALL ARRMC(3,1,DETAROH,-SCALR,DETAROH) 3969C WRITE(6,*) 'SCALR ',SCALR 3970C WRITE(6,*) 'DetaROH = ',DETAROH 3971 CALL ARRMC(3,1,ROH,1.,ROHOLD) 3972C SAVE THE PREVIOUS ROH 3973 CALL ARRAD(3,1,ROHOLD,DETAROH,ROH) 3974C ROH=ROHOLD+DETAROH 3975 CALL HUBERARC(ROH,A) 3976C NEW MATRIX 3977 CALL MATMULT(3,3,3,NATM,A,VEC1,VEC1P) 3978 CALL ARRPS(3,NATM,VEC1P,VEC2,DIS) 3979C A*V1 - V2 3980 RMS=dosq(3*NATM,DIS) 3981C 3982C WRITE(6,*) 'RMS NEW AND OLD', RMS,RMS1 3983 IF (RMS.LT.RMS1) THEN 3984 RMS1=RMS 3985 NREF=NREF+1 3986 GOTO 50 3987 ELSE IF(SCALR.GT.1E-4) THEN 3988 CALL ARRMC(3,1,DETAROH,-1./SCALR,DETAROH) 3989 SCALR = SCALR*0.5 3990 CALL ARRMC(3,1,ROHOLD,1.,ROH) 3991 CALL HUBERARC(ROH,A) 3992 GOTO 30 3993 END IF 3994 3995C WRITE(6,'(A)') '0' 3996C WRITE(6,*) 'That is the final rotation result.' 3997C WRITE(6,'(A)') '0' 3998 CALL ARRMC(3,1,ROHOLD,1.,ROH) 3999 CALL HUBERARC(ROH,A) 4000 4001 END 4002c 4003c 4004 SUBROUTINE REFRT ( NATM, X1, X2, A, T, VT, DIS ) 4005C DIMENSION X1(3,NATM), X2(3,NATM), A(3,3), T(3) 4006C DIMENSION VT(3*NATM,6),DIS(3,NATM) 4007C A subroutine to refine the superimpose matrix and vector from MOL1 to 4008c MOL2 i.e. XYZ2(3*NATM) = A(3*3) * XYZ1(3*NATM) + T(3) 4009C Where NATM represents the number of the atoms in each molecule 4010c which will be superimposed. 4011c XYZ1(3,NATM) represents the coordinates of MOL1 4012c XYZ2(3,NATM) represents the coordinates of MOL2 4013C A(3*3) represents the input initial matrix and 4014c output matrix. 4015c T(3*3) represents the input initial translation vector 4016c output vector. 4017C VT(3*6*NATM) is a 6*NATM-dimensional buffer array 4018C DIS(3*NATM) is a 3*NATM-dimensional buffer array 4019c 4020C COMMON/SHIFS/ SCALR,SCALT 4021c DATA SCALR/1./,SCALT/1./ 4022C Note: There should be a scale of the shifts. i.e. Shift= scale * shifts 4023c Here SCALR is the rotation shiftscale 4024c Here SCALT is the translation shiftscale 4025c The initial and suggested SCALR is 60. 4026c The initial and suggested SCALt is 1. 4027C 4028C COMMON/RMS/ RMS,SMEAN,NREF,NREF1 4029C RMS is the final R.M.S between two molecules. 4030C SMEAN is the mean distance. 4031c NREF is the number of cycles of refinement 4032c NREF1 is not useful. 4033c They could be used to judge the SCALR and SCALS 4034c 4035c The program use translation linear to least square refinement. Testing 4036c shows that it works quite well. The function of this routine is similar 4037c to SUPOS but the r.m.s. is less and the orientation might be more dangerous. 4038c 4039c 4040c by Guoguang Lu 4041C 27/01/1989 4042C 4043 COMMON/RMS/ RMS,SMEAN,NREF1,NREF 4044c COMMON/SHIFS/ SCALR,SCALT 4045c DATA SCALR/1./,SCALT/1./,IAT/0/ 4046C 4047 DIMENSION X1(3,NATM), X2(3,NATM), A(3,3), T(3) 4048 DIMENSION VT(NATM*3,6),DIS(3,NATM) 4049 DIMENSION VVV(6,6), VV1(6) 4050C 4051 DIMENSION DA(3,3,3),DROH(3),DT(3),DETA(6) 4052 DIMENSION B1(3),B2(3),ROHOLD(3),ROH(3),TOLD(3) 4053 REAL EMT(3,3) 4054 EQUIVALENCE (DETA,DROH) 4055 EQUIVALENCE (DETA(4),DT) 4056 data emt /1.,0.,0.,0.,1.,0.,0.,0.,1./ 4057C 4058 DEG=180./3.14159265359 4059 SCALR = 1. 4060 SCALT = 1. 4061 iat = 0. 4062 NREF=0 4063 epsilon = 1. 4064C 4065 CALL MTOHUBER(A,ROH,B1) 4066c 4067 CALL RTMOV(NATM,X1,A,T,DIS) 4068 CALL ARRPS(3,NATM,DIS,X2,DIS) 4069 RMS = DOSQ(NATM*3, DIS) 4070 RMS1 = RMS 4071 DO 5 I = 1, 3 4072 DO 5 J = 1, NATM 4073 5 CALL ARRGIVE(3,EMT(1,I),VT((J-1)*3+1,I+3)) 4074C 4075C Refine th1, th2, th3, tx, ty and tz 4076c 407710 CONTINUE 4078c 4079 if ( (nref.gt.99) .or. (abs(epsilon).le.1.0e-5) ) goto 25 4080C WRITE(6,*) 4081c WRITE(6,*) 'Cycle of refinement:',NREF 4082C WRITE(6,*) 'Matrix' 4083C WRITE(6,20) A 4084C20 FORMAT(3F10.5) 4085C WRITE(6,*) 'Eulerian Angle:',ROH 4086C WRITE(6,*) 'Translation: ',T 4087C WRITE(6,*) 'SIGMA (distance^2)=',RMS 4088C 4089C compute the normal matrix 4090 DO I =1, 3 4091 CALL DRVROHTH(I,ROH,DA(1,1,I)) 4092 CALL MATMULT(3,3,3,NATM,DA(1,1,I),X1,VT(1,I)) 4093C 4094c CALL ARRMC(3,1,B1,0.,B1) 4095c B1(I) = 1. 4096c DO J = 1, NATM 4097c CALL ARRMC(3,1,B1, 1., VT((J-1)*3+1,I+3) ) 4098c END DO 4099C 4100 END DO 4101C 4102 CALL LSQEQ(3*NATM,6,VT,DIS,DETA,VVV,VV1) 4103C 4104C SHIFT = SHIFT * SCALE 4105 CALL ARRMC(3,1,DROH,-SCALR*deg,DROH) 4106 CALL ARRMC(3,1,DT,-SCALT,DT) 4107C TYPE *, 'Shift ROH',DROH 4108C TYPE *, 'Shift T ',DT 4109C SAVE THE OLD ANGLE AND TRANSLATION VECTOR 4110 CALL ARRMC(3,1,ROH,1.,ROHOLD) 4111 CALL ARRMC(3,1,T,1.,TOLD) 4112C 4113 CALL ARRAD(3,1,ROH,DROH,ROH) 4114 CALL ARRAD(3,1,T,DT,T) 4115C 4116 CALL HUBER(ROH,A) 4117C 4118 CALL RTMOV(NATM,X1,A,T,DIS) 4119 CALL ARRPS(3,NATM,DIS,X2,DIS) 4120C 4121 RMS = DOSQ(NATM*3, DIS) 4122c TYPE *, 'rms new and old ',rms, rms1 4123C 4124 IF (RMS.LT.RMS1) THEN 4125 epsilon = (rms1-rms)/rms1 4126 NREF = NREF+1 4127 RMS1=RMS 4128 GOTO 10 4129 END IF 4130C 4131 25 continue 4132 CALL ARRMC(3,1,ROHOLD,1.,ROH) 4133 CALL ARRMC(3,1,TOLD,1.,T) 4134 CALL HUBER(ROH,A) 4135C 4136 CALL RTMOV(NATM,X1,A,T,DIS) 4137 CALL ARRPS(3,NATM,DIS,X2,DIS) 4138C 4139 ATM = NATM 4140 SMEAN = 0 4141 DO I=1, NATM 4142 D = VEM(3,DIS(1,I)) 4143 SMEAN = SMEAN + D 4144 END DO 4145 SMEAN = SMEAN /ATM 4146C 4147c SMEAN = SIGMA (DISTANCE) / N 4148C i i 4149 RMS = SQRT( RMS1/ATM ) 4150C 4151c RMS = SQRT( SIGMA (DISTANCE^2) / N ) 4152C i i 4153c 4154c 4155 END 4156c 4157c 4158 SUBROUTINE REFRTFIN ( NATM, X1, X2, A, T, VT, DIS ) 4159C DIMENSION X1(3,NATM), X2(3,NATM), A(3,3), T(3) 4160C DIMENSION VT(3*NATM,6),DIS(3,NATM) 4161C A subroutine to refine the superimpose matrix and vector from MOL1 to 4162c MOL2 i.e. XYZ2(3*NATM) = A(3*3) * XYZ1(3*NATM) + T(3) 4163C Where NATM represents the number of the atoms in each molecule 4164c which will be superimposed. 4165c XYZ1(3,NATM) represents the coordinates of MOL1 4166c XYZ2(3,NATM) represents the coordinates of MOL2 4167C A(3*3) represents the input initial matrix and 4168c output matrix. 4169c T(3*3) represents the input initial translation vector 4170c output vector. 4171C VT(3*6*NATM) is a 6*NATM-dimensional buffer array 4172C DIS(3*NATM) is a 3*NATM-dimensional buffer array 4173c 4174C COMMON/SHIFS/ SCALR,SCALT 4175c DATA SCALR/1./,SCALT/1./ 4176C Note: There should be a scale of the shifts. i.e. Shift= scale * shifts 4177c Here SCALR is the rotation shiftscale 4178c Here SCALT is the translation shiftscale 4179c The initial and suggested SCALR is 60. 4180c The initial and suggested SCALt is 1. 4181C 4182C COMMON/RMS/ RMS,SMEAN,NREF,NREF1 4183C RMS is the final R.M.S between two molecules. 4184C SMEAN is the mean distance. 4185c NREF is the number of cycles of refinement 4186c NREF1 is not useful. 4187c They could be used to judge the SCALR and SCALS 4188c 4189c The routine uses translation linear to least square refinement. Testing 4190c shows that it works quite well. The function of this routine is similar 4191c to SUPOS but the r.m.s. is less and the orientation might be more dangerous. 4192c 4193c 4194c by Guoguang Lu 4195C 27/01/1989 4196C 4197 COMMON/RMS/ RMS,SMEAN,NREF1,NREF 4198c COMMON/SHIFS/ SCALR,SCALT 4199C 4200 DIMENSION X1(3,NATM), X2(3,NATM), A(3,3), T(3) 4201 DIMENSION VT(NATM*3,6),DIS(3,NATM) 4202 DIMENSION VVV(6,6), VV1(6) 4203C 4204 DIMENSION DA(3,3,3),DROH(3),DT(3),DETA(6) 4205 DIMENSION B1(3),B2(3),ROHOLD(3),ROH(3),TOLD(3) 4206 EQUIVALENCE (DETA,DROH) 4207 EQUIVALENCE (DETA(4),DT) 4208c... data statements. Separate declaration and init required for f2c 4209 DATA SCALR/1./,SCALT/1./,IAT/0/ 4210C 4211 DEG=180./3.14159265359 4212 NREF=0 4213C 4214 CALL MTOHUBERARC(A,ROH,B1) 4215c 4216 CALL RTMOV(NATM,X1,A,T,DIS) 4217 CALL ARRPS(3,NATM,DIS,X2,DIS) 4218 RMS = DOSQ(NATM*3, DIS) 4219 RMS1 = RMS 4220C 4221C Refine th1, th2, th3, tx, ty and tz 4222c 422310 CONTINUE 4224c 4225C WRITE(6,*) 4226c WRITE(6,*) 'Cycle of refinement:',NREF 4227C WRITE(6,*) 'Matrix' 4228C WRITE(6,20) A 4229C20 FORMAT(3F10.5) 4230C WRITE(6,*) 'Eulerian Angle:',ROH 4231C WRITE(6,*) 'Translation: ',T 4232C WRITE(6,*) 'SIGMA (distance^2)=',RMS 4233C 4234C compute the normal matrix 4235 DO I =1, 3 4236 CALL DRVROHTHARC(I,ROH,DA(1,1,I)) 4237 CALL MATMULT(3,3,3,NATM,DA(1,1,I),X1,VT(1,I)) 4238C 4239 CALL ARRMC(3,1,B1,0.,B1) 4240 B1(I) = 1. 4241 DO J = 1, NATM 4242 CALL ARRMC(3,1,B1, 1., VT((J-1)*3+1,I+3) ) 4243 END DO 4244C 4245 END DO 4246C 4247 CALL LSQEQ(3*NATM,6,VT,DIS,DETA,VVV,VV1) 4248C 4249C SHIFT = SHIFT * SCALE 4250C CALL ARRMC(3,1,DROH,-SCALR*deg,DROH) 425130 continue 4252 CALL ARRMC(3,1,DROH,-SCALR,DROH) 4253 CALL ARRMC(3,1,DT,-SCALT,DT) 4254C TYPE *, 'scalr,scalt',scalr,scalt 4255C TYPE *, 'Shift ROH',DROH 4256C TYPE *, 'Shift T ',DT 4257C SAVE THE OLD ANGLE AND TRANSLATION VECTOR 4258 CALL ARRMC(3,1,ROH,1.,ROHOLD) 4259 CALL ARRMC(3,1,T,1.,TOLD) 4260C 4261 CALL ARRAD(3,1,ROH,DROH,ROH) 4262 CALL ARRAD(3,1,T,DT,T) 4263C 4264 CALL HUBERARC(ROH,A) 4265C 4266 CALL RTMOV(NATM,X1,A,T,DIS) 4267 CALL ARRPS(3,NATM,DIS,X2,DIS) 4268C 4269 RMS = DOSQ(NATM*3, DIS) 4270C TYPE *, 'refrt rms new and old ',rms, rms1 4271C 4272 IF (RMS.LT.RMS1) THEN 4273 NREF = NREF+1 4274 RMS1=RMS 4275 GOTO 10 4276 else if (scalr.gt.1.e-4) then 4277 CALL ARRMC(3,1,DROH,-1./SCALR,DROH) 4278 CALL ARRMC(3,1,DT,-1./SCALT,DT) 4279 scalr = scalr*0.5 4280 scalt = scalt*0.5 4281 CALL ARRMC(3,1,ROHOLD,1.,ROH) 4282 CALL ARRMC(3,1,TOLD,1.,T) 4283 CALL HUBERARC(ROH,A) 4284 goto 30 4285c 4286 END IF 4287C 4288 CALL ARRMC(3,1,ROHOLD,1.,ROH) 4289 CALL ARRMC(3,1,TOLD,1.,T) 4290 CALL HUBERARC(ROH,A) 4291C 4292 CALL RTMOV(NATM,X1,A,T,DIS) 4293 CALL ARRPS(3,NATM,DIS,X2,DIS) 4294C 4295 ATM = NATM 4296 SMEAN = 0 4297 DO I=1, NATM 4298 D = VEM(3,DIS(1,I)) 4299 SMEAN = SMEAN + D 4300 END DO 4301 SMEAN = SMEAN /ATM 4302C 4303c SMEAN = SIGMA (DISTANCE) / N 4304C i i 4305 RMS = SQRT( RMS1/ATM ) 4306C 4307c RMS = SQRT( SIGMA (DISTANCE^2) / N ) 4308C i i 4309c 4310c 4311 END 4312c 4313c 4314 SUBROUTINE REFRTFIN1 ( NATM, X1, X2, A, T, VT, DIS, DIST ) 4315C DIMENSION X1(3,NATM), X2(3,NATM), A(3,3), T(3) 4316C DIMENSION VT(NATM,6),DIS(3,NATM) 4317C A subroutine to refine the superimpose matrix and vector from MOL1 to 4318c MOL2 i.e. XYZ2(3*NATM) = A(3*3) * XYZ1(3*NATM) + T(3) 4319C Where NATM represents the number of the atoms in each molecule 4320c which will be superimposed. 4321c XYZ1(3,NATM) represents the coordinates of MOL1 4322c XYZ2(3,NATM) represents the coordinates of MOL2 4323C A(3*3) represents the input initial matrix and 4324c output matrix. 4325c T(3*3) represents the input initial translation vector 4326c output vector. 4327C VT(6*NATM) is a 6*NATM-dimensional buffer array 4328C DIS(3*NATM) is a 3*NATM-dimensional buffer array 4329C DIST(NATM) is a 3*NATM-dimensional buffer array 4330c 4331C COMMON/SHIFS/ SCALR,SCALT 4332c DATA SCALR/1./,SCALT/1./ 4333C Note: There should be a scale of the shifts. i.e. Shift= scale * shifts 4334c Here SCALR is the rotation shiftscale 4335c Here SCALT is the translation shiftscale 4336c The initial and suggested SCALR is 60. 4337c The initial and suggested SCALt is 1. 4338C 4339C COMMON/RMS/ RMS,SMEAN,NREF,NREF1 4340C RMS is the final R.M.S between two molecules. 4341C SMEAN is the mean distance. 4342c NREF is the number of cycles of refinement 4343c NREF1 is not useful. 4344c They could be used to judge the SCALR and SCALS 4345c 4346c The routine uses translation linear to least square refinement. Testing 4347c shows that it works quite well. The function of this routine is similar 4348c to SUPOS but the r.m.s. is less and the orientation might be more dangerous. 4349c 4350c 4351c by Guoguang Lu 4352C 27/01/1989 4353C 4354 COMMON/RMS/ RMS,SMEAN,NREF1,NREF 4355c COMMON/SHIFS/ SCALR,SCALT 4356c DATA SCALR/1./,SCALT/1./,IAT/0/ 4357C 4358 DIMENSION X1(3,NATM), X2(3,NATM), A(3,3), T(3) 4359 DIMENSION VT(NATM,6),DIS(3,NATM),DIST(NATM) 4360 DIMENSION VVV(6,6), VV1(6) 4361C 4362 DIMENSION DA(3,3,3),DROH(3),DT(3),DETA(6) 4363 DIMENSION B1(3),B2(3),ROHOLD(3),ROH(3),TOLD(3) 4364 EQUIVALENCE (DETA,DROH) 4365 EQUIVALENCE (DETA(4),DT) 4366C 4367c DEG=180./3.14159265359 4368 scalr = 1.0 4369 scalt = 1.0 4370 NREF=0 4371C 4372 CALL MTOHUBERarc(A,ROH,B1) 4373c 4374 CALL RTMOV(NATM,X1,A,T,DIS) 4375 CALL ARRPS(3,NATM,DIS,X2,DIS) 4376 DO I =1, NATM 4377 DIST(I) = VEM(3,DIS(1,I)) 4378 END DO 4379 RMS = DOSQ(NATM*3, DIS) 4380 RMS1 = RMS 4381C 4382C Refine th1, th2, th3, tx, ty and tz 4383c 438410 CONTINUE 4385c 4386c WRITE(6,*) 4387c WRITE(6,*) 'Cycle of refinement:',NREF 4388c WRITE(6,*) 'Matrix' 4389c WRITE(6,20) A 4390c20 FORMAT(3F10.5) 4391c WRITE(6,*) 'Eulerian Angle:',ROH 4392c WRITE(6,*) 'Translation: ',T 4393c WRITE(6,*) 'SIGMA (distance^2)=',RMS 4394C 4395C compute the normal matrix 4396 DO I =1, 3 4397 CALL DRVROHTHarc(I,ROH,DA(1,1,I)) 4398C 4399 DO J = 1, NATM 4400 CALL MATMULT(3,3,3,1,DA(1,1,I),X1(1,J),B1) 4401 VT(J,I) = POIMULT(3,3,DIS(1,J),B1) / DIST(J) 4402 VT(J,I+3) = DIS(I,J) / DIST(J) 4403 END DO 4404C 4405 END DO 4406 CALL LSQEQ(NATM,6,VT,DIST,DETA,VVV,VV1) 4407C TYPE *, 'Shift ', DETA 4408C SHIFT = SHIFT * SCALe 4409c CALL ARRMC(3,1,DROH,-SCALR*deg,DROH) 441030 continue 4411C TYPE *, 'SCALR, SCALT ', SCALR, SCALT 4412 CALL ARRMC(3,1,DROH,-SCALR,DROH) 4413 CALL ARRMC(3,1,DT,-SCALT,DT) 4414C TYPE *, 'Shift ROH',DROH 4415C TYPE *, 'Shift T ',DT 4416C SAVE THE OLD ANGLE AND TRANSLATION VECTOR 4417 CALL ARRMC(3,1,ROH,1.,ROHOLD) 4418 CALL ARRMC(3,1,T,1.,TOLD) 4419C 4420 CALL ARRAD(3,1,ROH,DROH,ROH) 4421 CALL ARRAD(3,1,T,DT,T) 4422C 4423 CALL HUBERarc(ROH,A) 4424C 4425 CALL RTMOV(NATM,X1,A,T,DIS) 4426 CALL ARRPS(3,NATM,DIS,X2,DIS) 4427 DO I =1, NATM 4428 DIST(I) = VEM(3,DIS(1,I)) 4429 END DO 4430C 4431 RMS = DOSQ(NATM*3, DIS) 4432C TYPE *, 'rms new and old ',rms, rms1 4433C 4434 IF (RMS.LT.RMS1) THEN 4435 NREF = NREF+1 4436 RMS1=RMS 4437 GOTO 10 4438c else IF (RMS.ge.0) THEN 4439 else IF (scalr.Gt.1.e-4) THEN 4440 CALL ARRMC(3,1,DROH,-1./SCALR,DROH) 4441 CALL ARRMC(3,1,DT,-1./SCALT,DT) 4442 scalr = scalr*0.5 4443 scalt = scalt*0.5 4444C 4445 CALL ARRMC(3,1,ROHOLD,1.,ROH) 4446 CALL ARRMC(3,1,TOLD,1.,T) 4447 CALL HUBERarc(ROH,A) 4448 goto 30 4449 END IF 4450C 4451 CALL ARRMC(3,1,ROHOLD,1.,ROH) 4452 CALL ARRMC(3,1,TOLD,1.,T) 4453 CALL HUBERarc(ROH,A) 4454 CALL RTMOV(NATM,X1,A,T,DIS) 4455 CALL ARRPS(3,NATM,DIS,X2,DIS) 4456 DO I =1, NATM 4457 DIST(I) = VEM(3,DIS(1,I)) 4458 END DO 4459 4460c RMS = DOSQ(NATM*3, DIS) 4461c TYPE *, 'rms repeat and old ',rms, rms1 4462 4463 4464c TYPE *, 'ending rms ',rms 4465 ATM = NATM 4466 SMEAN = 0 4467 DO I=1, NATM 4468 D = VEM(3,DIS(1,I)) 4469 SMEAN = SMEAN + D 4470 END DO 4471 SMEAN = SMEAN /ATM 4472c 4473c SMEAN = SIGMA (DISTANCE) / N 4474C i i 4475 RMS = SQRT( RMS1/ATM ) 4476C 4477c RMS = SQRT( SIGMA (DISTANCE^2) / N ) 4478C i i 4479c 4480c 4481 END 4482c 4483c 4484 function gg_res3to1(rin) 4485 character gg_res3to1*1,rin*3 4486 character*3 res(20),abr(20) 4487c... data statements. Separate declaration and init required for f2c 4488 data RES /'ALA','GLU','GLN','ASP','ASN','LEU','GLY' 4489 1 ,'LYS','SER','VAL','ARG','THR' 4490 2 ,'PRO','ILE','MET','PHE','TYR','CYS','TRP','HIS'/ 4491 data ABR /'A ','E ','Q ','D ','N ','L ','G ' 4492 1 ,'K ','S ','V ','R ','T ' 4493 2 ,'P ','I ','M ','F ','Y ','C ','W ','H '/ 4494c 4495 do i = 1, 20 4496 if (rin.eq.res(i)) then 4497 gg_res3to1 = abr(i)(1:1) 4498 return 4499 end if 4500 end do 4501 if (rin.eq.' ') then 4502 gg_res3to1 = ' ' 4503 return 4504 end if 4505c write(6,*) 'Warning: no such residues ',rin 4506 gg_res3to1 = ' ' 4507 end 4508c 4509c 4510 subroutine rotvsvec(vec,vkapa,a) 4511c a subroutine to give a matrix when you know there is a rotation 4512c which angle is vkapa along the vector vec as axis. 4513c vec is the input 3-dimensional vector where the rotation axis is 4514c vkapa is the rotating angle 4515c a is the output 3*3 rotating matrix. 4516c 4517 real vec(3),a(3,3) 4518 external sind 4519c WRITE(6,*) 'vec,vkapa' 4520C WRITE(6,*) vec,vkapa 4521c There is a strange bug on Alpha, ---- cosd(vkapa) has the wrong sign!!!! 4522c Instead I have to have cos(vkapa*pai/180.) 4523 pai = 3.141592654 4524c 4525 sinkapa = sind(vkapa) 4526 coskapa = cos(vkapa*pai/180.) 4527c 4528 a(1,1) = (1.-vec(1)*vec(1))*coskapa+vec(1)*vec(1) 4529 a(2,2) = (1.-vec(2)*vec(2))*coskapa+vec(2)*vec(2) 4530 a(3,3) = (1.-vec(3)*vec(3))*coskapa+vec(3)*vec(3) 4531c 4532 a(1,2) = vec(1)*vec(2)*(1.-coskapa) -vec(3)*sinkapa 4533 a(2,1) = vec(1)*vec(2)*(1.-coskapa) +vec(3)*sinkapa 4534c 4535 a(1,3) = vec(3)*vec(1)*(1.-coskapa) +vec(2)*sinkapa 4536 a(3,1) = vec(3)*vec(1)*(1.-coskapa) -vec(2)*sinkapa 4537c 4538 a(2,3) = vec(2)*vec(3)*(1.-coskapa) -vec(1)*sinkapa 4539 a(3,2) = vec(2)*vec(3)*(1.-coskapa) +vec(1)*sinkapa 4540c 4541 end 4542c 4543c 4544 SUBROUTINE RTMOV(NATM,XIN,A,T,XOUT) 4545C A routine to compute: XOUT = A * XIN + T 4546C where XIN is the input 3*NATM array 4547C XOUT is the Output 3*NATM array 4548c A(3,3) is a 3*3 rotation matrix 4549c T(3) is a 3-dimensional translation vector. 4550c DIMENSION XIN(3,NATM),XOUT(3,NATM),A(3,3),T(3) 4551C 4552 DIMENSION XIN(3,NATM),XOUT(3,NATM),A(3,3),T(3) 4553 CALL MATMULT(3,3,3,NATM,A,XIN,XOUT) 4554 DO I = 1, NATM 4555 CALL ARRAD(3,1,XOUT(1,I),T,XOUT(1,I)) 4556 END DO 4557 END 4558c 4559c 4560 function screw(a,t) 4561c dimension a(3,3),t(3),POL(3),xyz0(3) 4562C this is a routine to calculate how much the movement is along the 4563c the rotation axis. 4564c When there is a non-crystallographic rotation and translation 4565c x2 = A*x1 + t 4566c where A(3,3) is the 3*3 matrix 4567c t(3) is translation vector. 4568c the routine output POL(3) (psi, phi and kappa), go and xyz0(3) 4569c where psi (POL(1)) is the angle between the rotation axis and y axis 4570c phi is the angle between X axis and the image of the rotation axis 4571c kappa is the rotating angle. 4572c go is the screw amount along the axis direction 4573c Guoguang 900927 4574c 4575 4576 dimension a(3,3),t(3) 4577 dimension vec(3) 4578c 4579 call mtovec(a,vec,vkapa) 4580 if (vkapa.lt.1e-6) then 4581 screw = vem(3,t) 4582 end if 4583c 4584 if (vem(3,t).lt.(1.e-6)) then 4585 screw = 0. 4586 return 4587 end if 4588c vec is the direction of the axis vector. 4589 screw = poimult(3,3,t,vec) 4590 end 4591c 4592c 4593 subroutine sectime(sec,nhour,nmin,sres) 4594c input seconds and output the hour and minute and residue second 4595 nmin = 0 4596 nhour = 0 4597 if (sec.gt.60.) then 4598 smin = sec/60. 4599 nmin = int(smin) 4600 sres = sec - nmin*60 4601 end if 4602 if (nmin.gt.60) then 4603 nmin0 = nmin 4604 nhour = int(float(nmin)/60.) 4605 nmin = nmin0 - nhour*60 4606 end if 4607 end 4608c 4609c 4610 subroutine seekaxis(a,t,cen1,cen2,POL,go,xyz0) 4611c dimension a(3,3),t(3),POL(3),xyz0(3) 4612C this is a routine to calculate how much the movement is along the 4613c the rotation axis and where is the rotation axis from a rigid movement 4614c matrix. 4615c When there is a non-crystallographic rotation and translation 4616c x2 = A*x1 + t 4617c where A(3,3) is the 3*3 matrix 4618c t(3) is translation vector. 4619c cen(3), center of the two mass molecule 4620c the routine output POL(3) (psi, phi and kappa), go and xyz0(3) 4621c where psi (POL(1)) is the angle between the rotation axis and y axis 4622c phi is the angle between X axis and the image of the rotation axis 4623c kappa is the rotating angle. 4624c go is the screw amount along the axis direction 4625c Guoguang 900927 4626c xyz0(3) is one point where the axis pass by. This point should be 4627c a project point of CEN on the axis line, where CEN is center of 4628c two input molecules (CEN1+CEN2)/2. 4629c 4630 4631 dimension a(3,3),t(3),POL(3),xyz0(3),cen(3),cen1(3),cen2(3) 4632 dimension b1(3),b2(3),b3(3) 4633 dimension rot(3,3), rot1(3,3) 4634 dimension b4(3),b5(3),cenp1(3),cenp2(3),xyzp(3) 4635 dimension vec(3) 4636 logical*1 zax 4637 external tand 4638c 4639 call arrmc(3,1,t,0.,xyz0) 4640 call mtopolorz(a,POL,b1) 4641 if (abs(POL(3)).lt.1e-6) then 4642 go = vem(3,t) 4643 end if 4644c 4645 if (vem(3,t).lt.(1.e-6)) then 4646 go = 0. 4647 return 4648 end if 4649c 4650 call arrad(3,1,cen1,cen2,b4) 4651 call arrmc(3,1,b4,0.5,cen) 4652c 4653c b1(2) = cosd(POL(1)) 4654c b1(1) = sind(POL(1))*cosd(POL(2)) 4655c b1(3) = - sind(POL(1))*sind(POL(2)) 4656 call mtovec(a,vec,vkapa) 4657 call arrgive(3,vec,b3) 4658c b1 and vec is the direction of the axis vector. 4659 go = poimult(3,3,t,b3) 4660c now define a new coordinate system 4661c X0 = rot * X' 4662c where Z' is the rotating axis. 4663c and X'is in the Z'/cen2-cen1 plane 4664c Y' is perpendicular to X' and Z' 4665c so b1 is Z' vector. 4666 call arrps(3,1,cen2,cen1,b4) 4667 call veccrsmlt(b3,b4,b5) 4668 temp = vem(3,b5) 4669 if (temp.gt.(1.e-6)) then 4670 zax = .false. 4671 call arrmc(3,1,b5,1./temp,b2) 4672c here b2 is X' vector 4673 call veccrsmlt(b2,b3,b1) 4674c b3 is Y' vector 4675c 4676 call arrgive(3,b2,rot(1,2)) 4677 call arrgive(3,b3,rot(1,3)) 4678 call arrgive(3,b1,rot(1,1)) 4679 else 4680 call elize(3,rot) 4681 zax = .true. 4682 end if 4683c 4684 call antiarr(3,3,rot,rot1) 4685c rot1 = rot**-1 4686c 4687c b1 is the translation in new coordinate system 468811 call matmult(3,3,3,1,rot1,t,b1) 4689 if (abs(b1(3)-go).gt.1e-6) then 4690 if (zax) then 4691 rot(3,3) = -1. 4692 rot(2,2) = -1. 4693 rot1(3,3) = -1. 4694 rot1(2,2) = -1. 4695 zax = .false. 4696 goto 11 4697 else 4698 write(6,*) 'go ', go 4699 write(6,*) 't ',t 4700 write(6,*) 'b1 ',b1 4701 write(6,*) 'rot1' 4702 write(6,'(3f12.5)') rot1 4703 stop 'guoguang, you are wrong to get go' 4704 end if 4705 end if 4706c 4707c take cen as origin. take cen1 cen2 into new coordinates as cenp1,cenp2 4708 call arrps(3,1,cen1,cen,b1) 4709 call matmult(3,3,3,1,rot1,b1,cenp1) 4710 call arrps(3,1,cen2,cen,b2) 4711 call matmult(3,3,3,1,rot1,b2,cenp2) 4712c 4713 if (cenp1(2).gt.1.e-4) print *, 'Warning cenp1',cenp1 4714 if (cenp2(2).gt.1.e-4) print *, 'Warning cenp2',cenp2 4715c the passing point in new system 4716 xyzp(1) = 0. 4717 xyzp(3) = 0. 4718 if (abs(180.-vkapa).lt.0.001) then 4719 xyzp(2) = 0. 4720 else 4721 xyzp(2) = -cenp1(1)/tand(vkapa/2.) 4722 end if 4723c the passing point in old system 4724 call matmult(3,3,3,1,rot,xyzp,b4) 4725 call arrad(3,1,b4,cen,xyz0) 4726c 4727 end 4728c 4729c 4730 subroutine seekscrew(a,t,POL,go,xyz0) 4731c dimension a(3,3),t(3),POL(3),xyz0(3) 4732C this is a routine to calculate how much the movement is along the 4733c the rotation axis. 4734c When there is a non-crystallographic rotation and translation 4735c x2 = A*x1 + t 4736c where A(3,3) is the 3*3 matrix 4737c t(3) is translation vector. 4738c the routine output POL(3) (psi, phi and kappa), go and xyz0(3) 4739c where psi (POL(1)) is the angle between the rotation axis and y axis 4740c phi is the angle between X axis and the image of the rotation axis 4741c kappa is the rotating angle. 4742c go is the screw amount along the axis direction 4743c Guoguang 900927 4744c xyz0(3) is one point where the axis pass by. 4745c 4746 4747 dimension a(3,3),t(3),POL(3),xyz0(3) 4748 dimension b1(3),b2(3),b3(3) 4749 dimension rot(3,3), rot1(3,3),rot3(3,3),rot2(3,3) 4750 dimension b4(3) 4751 dimension vec(3) 4752 logical*1 zax 4753 external acosd, cosd 4754c 4755 call arrmc(3,1,t,0.,xyz0) 4756 call mtopolors(a,POL,b1) 4757 if (abs(POL(3)).lt.1e-6) then 4758 go = vem(3,t) 4759 end if 4760c 4761 if (vem(3,t).lt.(1.e-6)) then 4762 go = 0. 4763 return 4764 end if 4765c 4766c b1(2) = cosd(POL(1)) 4767c b1(1) = sind(POL(1))*cosd(POL(2)) 4768c b1(3) = - sind(POL(1))*sind(POL(2)) 4769 call mtovec(a,b1,vkapa) 4770 call arrmc(3,1,b1,1.,vec) 4771c b1 and vec is the direction of the axis vector. 4772 go = poimult(3,3,t,b1) 4773c 4774 b2(1) = 0. 4775 b2(2) = 0. 4776 b2(3) = 1. 4777c now decide a new coordinate system 4778c X0 = rot * X' 4779c where Z' is the rotating axis. 4780c and X'is in the XY plane. 4781c so b1 is Z' vector. 4782 call veccrsmlt(b1,b2,b3) 4783 temp = vem(3,b3) 4784 if (temp.gt.(1.e-6)) then 4785 zax = .false. 4786 call arrmc(3,1,b3,1./temp,b2) 4787c here b2 is X' vector 4788 call veccrsmlt(b1,b2,b3) 4789c b3 is Y' vector 4790c 4791 call arrmc(3,1,b2,1.,rot(1,1)) 4792 call arrmc(3,1,b3,1.,rot(1,2)) 4793 call arrmc(3,1,b1,1.,rot(1,3)) 4794 else 4795 call elize(3,rot) 4796 zax = .true. 4797 end if 4798c 4799 call antiarr(3,3,rot,rot1) 4800c rot1 = rot**-1 4801c 4802c b1 is the translation in new coordinate system 480311 call matmult(3,3,3,1,rot1,t,b1) 4804 if (abs(b1(3)-go).gt.1e-6) then 4805 if (zax) then 4806 rot(3,3) = -1. 4807 rot(2,2) = -1. 4808 rot1(3,3) = -1. 4809 rot1(2,2) = -1. 4810 zax = .false. 4811 goto 11 4812 else 4813 write(6,*) 'go ', go 4814 write(6,*) 't ',t 4815 write(6,*) 'b1 ',b1 4816 write(6,*) 'rot1' 4817 write(6,'(3f12.5)') rot1 4818 stop 'guoguang, you are wrong to get go' 4819 end if 4820 end if 4821c 4822 if (vem(2,b1).lt.1e-6) then 4823c the axes pass by the origin. 4824 call arrmc(3,1,t,0.,xyz0) 4825 return 4826 end if 4827c 4828c set a new coordinate system to make axis x'' is the image of t in new 4829c x'y' plane. and X' = rot2 * X'' 4830c 4831c 4832 b1(3) = 0. 4833 call elize(3,rot2) 4834 call arrmc(3,1,b1, 1./vem(3,b1), rot2(1,1) ) 4835 call veccrsmlt(rot2(1,3),rot2(1,1),rot2(1,2)) 4836 call antiarr(3,3,rot2,rot3) 4837c x'' = rot3 * x' = rot3 * rot1 * x0 4838 call matmult(3,3,3,1,rot3,b1,b4) 4839 if ( abs(b4(1)-vem(2,b1)).gt.1e-3 .or. abs(b4(2)).gt.1e-3 .or 4840 1 . abs(b4(3)).gt.1e-3) THEN 4841 WRITE(6,*) 'b1',b1 4842 WRITE(6,*) 'b4',b4 4843 write(6,*) 'rot2' 4844 write(6,'(3f12.5)') rot2 4845 write(6,*) 'rot3' 4846 write(6,'(3f12.5)') rot3 4847 stop 'something wrong in rot2 or rot3' 4848 end if 4849 x0 = b4(1) 4850c 4851 b2(3) = 0. 4852 b2(1) = x0/2 4853 vkapa = POL(3) 4854 b2(2) = sqrt( (1+cosd(vkapa))/(1-cosd(vkapa)) ) * b2(1) 4855 if (vkapa.lt.0.) b2(2) = - b2(2) 4856c b2 is what you want in x'' system 4857c the mathematics is very simple, isn't it? 4858c 4859c check result 4860 call arrps(3,1,b2,b4,b3) 4861 ang = acosd( poimult(3,3,b2,b3)/(vem(3,b2)*vem(3,b3)) ) 4862c if (abs(temp).gt.1) WRITE(6,*) temp 4863c ang = asind(temp) 4864 if (abs(ang-abs(vkapa)).gt.1e-2) then 4865 WRITE(6,*) 'something wrong in seekscrew' 4866 WRITE(6,*) 'angle ',ang 4867 WRITE(6,*) 'b1 ',b1 4868 WRITE(6,*) 'b2 ',b2 4869 WRITE(6,*) 'b3 ',b3 4870 WRITE(6,*) 'b4 ',b4 4871 end if 4872c now x0 = rot * x' = rot * rot2 * x'' = rot1 * x'' 4873 call matmult(3,3,3,3,rot,rot2,rot1) 4874 call matmult(3,3,3,1,rot1,b2,xyz0) 4875c to check the result. 4876 ds1 = dstpl2(vec,xyz0,t,b1,b4) 4877 call arrmc(3,1,t,0.,b3) 4878 ds2 = dstpl2(vec,xyz0,b3,b2,b4) 4879 df = abs(ds2-ds1) 4880c 4881 if (ds1.gt.1e-4) then 4882 if (df/ds1.gt.1e-4) 4883 + stop 'distance between 0-axis and t to axis is not same' 4884 4885 else 4886 if (ds2.gt.2e-4) 4887 + stop 'distance between 0-axis and t to axis is not same' 4888 end if 4889 4890 call arrps(3,1,b1,b2,b3) 4891 go1 = vem(3,b3) 4892 if (abs(abs(go)-go1) .gt. 1e-4) then 4893 WRITE(6,*) 'problem in go' 4894 WRITE(6,*) 'go ', go 4895 WRITE(6,*) 'screw', go1 4896 end if 4897c go is the real screw value along axis 4898 end 4899c 4900c 4901 subroutine seekscrewz(a,t,POL,go,xyz0) 4902c dimension a(3,3),t(3),POL(3),xyz0(3) 4903C this is a routine to calculate how much the movement is along the 4904c the rotation axis. 4905c When there is a non-crystallographic rotation and translation 4906c x2 = A*x1 + t 4907c where A(3,3) is the 3*3 matrix 4908c t(3) is translation vector. 4909c the routine output POL(3) (psi, phi and kappa), go and xyz0(3) 4910c where psi (POL(1)) is the angle between the rotation axis and y axis 4911c phi is the angle between X axis and the image of the rotation axis 4912c kappa is the rotating angle. 4913c go is the screw amount along the axis direction 4914c Guoguang 900927 4915c xyz0(3) is one point where the axis pass by. 4916c 4917 4918 dimension a(3,3),t(3),POL(3),xyz0(3) 4919 dimension b1(3),b2(3),b3(3) 4920 dimension rot(3,3), rot1(3,3),rot3(3,3),rot2(3,3) 4921 dimension b4(3) 4922 dimension vec(3) 4923 logical*1 zax 4924 external acosd, cosd 4925c 4926 call arrmc(3,1,t,0.,xyz0) 4927 call mtopolorz(a,POL,b1) 4928 if (abs(POL(3)).lt.1e-6) then 4929 go = vem(3,t) 4930 end if 4931c 4932 if (vem(3,t).lt.(1.e-6)) then 4933 go = 0. 4934 return 4935 end if 4936c 4937c b1(2) = cosd(POL(1)) 4938c b1(1) = sind(POL(1))*cosd(POL(2)) 4939c b1(3) = - sind(POL(1))*sind(POL(2)) 4940 call mtovec(a,b1,vkapa) 4941 call arrmc(3,1,b1,1.,vec) 4942c b1 and vec is the direction of the axis vector. 4943 go = poimult(3,3,t,b1) 4944c 4945 b2(1) = 0. 4946 b2(2) = 0. 4947 b2(3) = 1. 4948c now decide a new coordinate system 4949c X0 = rot * X' 4950c where Z' is the rotating axis. 4951c and X'is in the XY plane. 4952c so b1 is Z' vector. 4953 call veccrsmlt(b1,b2,b3) 4954 temp = vem(3,b3) 4955 if (temp.gt.(1.e-6)) then 4956 zax = .false. 4957 call arrmc(3,1,b3,1./temp,b2) 4958c here b2 is X' vector 4959 call veccrsmlt(b1,b2,b3) 4960c b3 is Y' vector 4961c 4962 call arrmc(3,1,b2,1.,rot(1,1)) 4963 call arrmc(3,1,b3,1.,rot(1,2)) 4964 call arrmc(3,1,b1,1.,rot(1,3)) 4965 else 4966 call elize(3,rot) 4967 zax = .true. 4968 end if 4969c 4970 call antiarr(3,3,rot,rot1) 4971c rot1 = rot**-1 4972c 4973c b1 is the translation in new coordinate system 497411 call matmult(3,3,3,1,rot1,t,b1) 4975 if (abs(b1(3)-go).gt.1e-6) then 4976 if (zax) then 4977 rot(3,3) = -1. 4978 rot(2,2) = -1. 4979 rot1(3,3) = -1. 4980 rot1(2,2) = -1. 4981 zax = .false. 4982 goto 11 4983 else 4984 write(6,*) 'go ', go 4985 write(6,*) 't ',t 4986 write(6,*) 'b1 ',b1 4987 write(6,*) 'rot1' 4988 write(6,'(3f12.5)') rot1 4989 stop 'guoguang, you are wrong to get go' 4990 end if 4991 end if 4992c 4993 if (vem(2,b1).lt.1e-6) then 4994c the axes pass by the origin. 4995 call arrmc(3,1,t,0.,xyz0) 4996 return 4997 end if 4998c 4999c set a new coordinate system to make axis x'' is the image of t in new 5000c x'y' plane. and X' = rot2 * X'' 5001c 5002c 5003 b1(3) = 0. 5004 call elize(3,rot2) 5005 call arrmc(3,1,b1, 1./vem(3,b1), rot2(1,1) ) 5006 call veccrsmlt(rot2(1,3),rot2(1,1),rot2(1,2)) 5007 call antiarr(3,3,rot2,rot3) 5008c x'' = rot3 * x' = rot3 * rot1 * x0 5009 call matmult(3,3,3,1,rot3,b1,b4) 5010 if ( abs(b4(1)-vem(2,b1)).gt.1e-3 .or. abs(b4(2)).gt.1e-3 .or 5011 1 . abs(b4(3)).gt.1e-3) THEN 5012 WRITE(6,*) 'b1',b1 5013 WRITE(6,*) 'b4',b4 5014 write(6,*) 'rot2' 5015 write(6,'(3f12.5)') rot2 5016 write(6,*) 'rot3' 5017 write(6,'(3f12.5)') rot3 5018 stop 'something wrong in rot2 or rot3' 5019 end if 5020 x0 = b4(1) 5021c 5022 b2(3) = 0. 5023 b2(1) = x0/2 5024 vkapa = POL(3) 5025 b2(2) = sqrt( (1+cosd(vkapa))/(1-cosd(vkapa)) ) * b2(1) 5026 if (vkapa.lt.0.) b2(2) = - b2(2) 5027c b2 is what you want in x'' system 5028c the mathematics is very simple, isn't it? 5029c 5030c check result 5031 call arrps(3,1,b2,b4,b3) 5032 ang = acosd( poimult(3,3,b2,b3)/(vem(3,b2)*vem(3,b3)) ) 5033c if (abs(temp).gt.1) WRITE(6,*) temp 5034c ang = asind(temp) 5035 if (abs(ang-abs(vkapa)).gt.1e-2) then 5036 WRITE(6,*) 'something wrong in seekscrew' 5037 WRITE(6,*) 'angle ',ang 5038 WRITE(6,*) 'b1 ',b1 5039 WRITE(6,*) 'b2 ',b2 5040 WRITE(6,*) 'b3 ',b3 5041 WRITE(6,*) 'b4 ',b4 5042 end if 5043c now x0 = rot * x' = rot * rot2 * x'' = rot1 * x'' 5044 call matmult(3,3,3,3,rot,rot2,rot1) 5045 call matmult(3,3,3,1,rot1,b2,xyz0) 5046c to check the result. 5047 ds1 = dstpl2(vec,xyz0,t,b1,b4) 5048 call arrmc(3,1,t,0.,b3) 5049 ds2 = dstpl2(vec,xyz0,b3,b2,b4) 5050 df = abs(ds2-ds1) 5051c 5052 if (ds1.gt.1e-4) then 5053 if (df/ds1.gt.1e-4) 5054 + stop 'distance between 0-axis and t to axis is not same' 5055 else 5056 if (ds2.gt.2e-4) 5057 + stop 'distance between 0-axis and t to axis is not same' 5058 end if 5059 5060 call arrps(3,1,b1,b2,b3) 5061 go1 = vem(3,b3) 5062 if (abs(abs(go)-go1) .gt. 1e-4) then 5063 WRITE(6,*) 'problem in go' 5064 WRITE(6,*) 'go ', go 5065 WRITE(6,*) 'screw', go1 5066 end if 5067c go is the real screw value along axis 5068 end 5069c 5070c 5071 SUBROUTINE SQSTLZ(M,P,PS,U,V) 5072C a subroutine to calculate the orientation matrix and the two axes 5073c from a symmetric matrix of an ellipse. by LGG 5074c M(2,2) is the input symmetric matrix of the ellipse. 5075C P(2,2) is one of the two orientation matrices of the ellipse 5076C PS(2,2) is the inverse matrix of the matrix PS(2,2) 5077C U, V are the lengths of the two axes of the ellipse. 5078C 5079 DIMENSION M(2,2),P(2,2),PS(2,2) 5080 REAL M,P,A,B,C,U,V,X1,X2 508110 FORMAT(1X,4F10.3) 5082 A=1 5083 B=-(M(1,1)+M(2,2)) 5084 C=M(1,1)*M(2,2)-M(2,1)*M(1,2) 5085 CALL SQUROOT(A,B,C,X1,X2) 5086 U=SQRT(1./X1) 5087 V=SQRT(1./X2) 5088 P(2,1)=(X1-M(1,1))/M(1,2) 5089 P(2,2)=(X2-M(1,1))/M(1,2) 5090 S1=SQRT(P(2,1)*P(2,1)+1) 5091 S2=SQRT(P(2,2)*P(2,2)+1) 5092 P(2,1)=P(2,1)/S1 5093 P(2,2)=P(2,2)/S2 5094 P(1,1)=1./S1 5095 P(1,2)=1./S2 5096 PS(1,2)=P(2,1) 5097 PS(2,1)=P(1,2) 5098 PS(1,1)=P(1,1) 5099 PS(2,2)=P(2,2) 5100 RETURN 5101 END 5102c 5103c 5104 subroutine spacegp(nunit,file,latnam,nsym,nrot,sym,nsp) 5105c In this subroutine the user gives a space group name, the subroutine reads the 5106c ccp4 symmetry library, then it outputs the symmetry operations 5107c Input 5108c nunit --- unit number of library file 5109c file--- symmetry libary file name 5110c latnam --- *14 character of a spacegroup name 5111c Output 5112c nsymm --- Total number of symmetry operations in this spacegroup 5113c nrot ---- only rotation symmetric operation in this spacegroup 5114c symm ---- a 3*4*nsym symmmetric operation c matrix. 5115c nsp --- space group number 5116c in common symdump, if dump is true, dump all the information 5117c maxop is the maximum number of symmetry operations 5118c library format 5119c155 6 2 R32 5120c X,Y,Z * -Y,X-Y,Z * Y-X,-X,Z 5121c Y,X,-Z * -X,Y-X,-Z * X-Y,-Y,-Z 5122c X+1/3,Y+2/3,Z+2/3 * -Y+1/3,X-Y+2/3,Z+2/3 * Y-X+1/3,-X+2/3,Z+2/3 5123c Y+1/3,X+2/3,-Z+2/3 * -X+1/3,Y-X+2/3,-Z+2/3 * X-Y+1/3,-Y+2/3,-Z+2/3 5124c X+2/3,Y+1/3,Z+1/3 * -Y+2/3,X-Y+1/3,Z+1/3 * Y-X+2/3,-X+1/3,Z+1/3 5125c Y+2/3,X+1/3,-Z+1/3 * -X+2/3,Y-X+1/3,-Z+1/3 * X-Y+2/3,-Y+1/3,-Z+1/3 5126 parameter (maxop=96) 5127 character file*80,latnam*14 5128 real sym(3,4,maxop) 5129 integer*4 nunit,nsp,nsym,nrot 5130 logical*1 dump 5131c COMMON/SYMDUMP/ DUMP 5132c 5133 character*80 key,str*40 5134c... data statements. Separate declaration and init required for f2c 5135 data dump /.false./ 5136c 5137c WRITE(6,*) 'nsp',nsp 5138 nsym=0 5139 if (nunit.eq.0) nunit=25 5140 if (file(1:1).eq.' ') call spstrunct(file) 5141 if (file(1:1).eq.' ') file='SYMOP' 5142 call ccpdpn(nunit,file(1:lenstr(file)),'READONLY','F',0,0) 514310 read(nunit,'(a)',end=220) key 5144 im = lenstr(latnam) 5145 in = lenstr(key) 5146 is = index(key(1:in),latnam(1:im)) 5147 if (is.eq.0) goto 10 5148c call redstrin(12,i1-1,key,npar) 5149 read(key(1:12),*,err=10,end=10) isp,iall,irot 5150c read(key(12:80),'(a)') latnam 5151c if (isp.ne.nsp) goto 10 5152 write(6,*) 'Space Group >>> ',Latnam,isp 5153 do i = 1, iall 5154 read(nunit,'(a)') key 5155 il = lenstr(key) 5156 key(il+1:il+1)='*' 5157 iend = 0 5158 if (dump) write(6,'(1x,a)') key(1:il) 515920 ist = iend+1 5160 iend = index(key(ist:il+1),'*')+IST-1 5161 str = ' ' 5162 str(1:iend-ist) = key(ist:iend-1) 5163 NSYM=NSYM+1 5164 IF (MATSYM(SYM(1,1,nsym),STR,ICOL).EQ.0) GOTO 502 5165 WRITE(6,*) 'Error in symop after column ',ICOL 5166 write(6,*) key 5167 write(6,*) str 5168 stop 'Check you file SYMOP' 5169502 continue 5170 if (dump) then 5171 write(6,'(4f8.4)') ((sym(j1,j2,nsym),j2=1,4),j1=1,3) 5172 end if 5173 if (iend.lt.il+1) goto 20 5174 if (i.eq.irot) nrot = nsym 5175 end do 5176 write(6,45) nsym,nrot 517745 format(1x,'Symmetric operation ----',6x, 5178 1 'Total: ',i3,6x,'Rotation:',i3) 5179 close (unit=nunit) 5180 return 5181210 WRITE(6,*) 'Error opening the SYMOP file ',FILE(1:LENSTR(FILE)) 5182220 WRITE(6,*) 'Space group',latnam(1:im), 5183 + ' was not found in SYMOP file' 5184 nsym = 0 5185 nrot = 0 5186 latnam = ' ' 5187 STOP 'Check your SYMOP file' 5188 END 5189c 5190c 5191 SUBROUTINE SPSTRUNCT(STRING) 5192 CHARACTER*(*) STRING 5193 LENS = LEN(STRING) 5194 IL = LENSTR(STRING) 5195C 51965 CONTINUE 5197 ISP = INDEX(STRING(1:IL),' ') 5198 IF (ISP.EQ.0.OR.ISP.GE.IL) RETURN 5199 STRING(ISP:IL-1) = STRING(ISP+1:IL) 5200 STRING(IL:IL) = ' ' 5201 IL = IL - 1 5202 GOTO 5 5203 END 5204c 5205c 5206 SUBROUTINE SQUROOT(A,B,C,X1,X2) 5207c Finds roots of a quadratic equation 5208 REAL A,B,C,X1,X2 5209 X1=(-B+SQRT(B*B-4*A*C))/(2*A) 5210 X2=(-B-SQRT(B*B-4*A*C))/(2*A) 5211 RETURN 5212 END 5213c 5214c 5215C a subroutine to calculate the statistics of a group of numbers 5216c subroutine statis(n,value,vmean,rms,dmean) 5217c dimension value(n) 5218c n is the number of value in the group 5219c value is the n-dimensional array whose statistics are required. 5220c vmean is the output mean value of the array value(n) 5221c rms is the output r.m.s in value(n) 5222c dmean is the output mean deviation of value(n) 5223 subroutine statis(n,value,vmean,rms,dmean) 5224 dimension value(n) 5225 vmean = 0 5226 do i =1, n 5227 vmean = vmean + value(i) 5228 end do 5229 vmean = vmean/n 5230c 5231 rms = 0. 5232 dmean = 0. 5233 do i =1, n 5234 temp = value(i) - vmean 5235 dmean = dmean + abs(temp) 5236 rms = rms + temp*temp 5237 end do 5238 rms = sqrt(rms/n) 5239 dmean = dmean / n 5240 end 5241c 5242c 5243 SUBROUTINE LGG_STRING(NUNIT,NCHA,TXT,NPAR,cont) 5244C A subroutine to write a character string to unit NUNIT using 5245c the spaces in the string to split the string into a list of 5246c parameters. 5247c NUNIT is the unit number. 5248c NCHA is the length of the character string. 5249c TXT is the character string. 5250C NPAR is the number of parameters in this string. 5251c cont is a flag which represent if this is a continue string or start 5252c .true. = continue 5253c .false. = start 5254c 5255c 5256 CHARACTER*(*) TXT 5257 logical*1 cont 5258 5259 if (.not.cont) then 5260 NPAR = 0 5261 REWIND (NUNIT) 5262 end if 5263 5264 jcha = 0 5265 if (ncha.gt.0) then 5266 if (txt(ncha:ncha).eq.'-') then 5267 jcha = 1 5268 txt(ncha:ncha)=' ' 5269 end if 5270 end if 5271 5272 I = 1 5273 10 IF (I.LE.NCHA-jcha) then 5274 if (TXT(I:I).NE.' ') THEN 5275 J = I 5276 20 CONTINUE 5277 IF ((J+1).GT.NCHA) THEN 5278 WRITE(NUNIT,'(A)') TXT(I:J) 5279 NPAR = NPAR + 1 5280 else if (TXT(J+1:J+1).EQ.' ') then 5281 WRITE(NUNIT,'(A)') TXT(I:J) 5282 NPAR = NPAR + 1 5283 ELSE 5284 J = J + 1 5285 GOTO 20 5286 END IF 5287 I = J + 1 5288 GOTO 10 5289 ELSE IF (I.LT.NCHA-jcha) THEN 5290 I = I + 1 5291 GOTO 10 5292 END IF 5293 endif 5294 5295 if (jcha.eq.1) TXT(NCHA:NCHA)='-' 5296C 5297c if (jcha.eq.0) REWIND(NUNIT) 5298C 5299 END 5300c 5301c 5302 SUBROUTINE SUPIM ( NATM, XYZ1, XYZ2, A, T ) 5303C DIMENSION XYZ1(3,NATM), XYZ2(3,NATM), A(3,3), T(3) 5304C A subroutine to give the superimpose matrix and vector from MOL1 to 5305c MOL2 i.e. XYZ2(3*NATM) = A(3*3) * XYZ1(3*NATM) + T(3) 5306C Where NATM represents the number of the atoms in each molecule 5307c which will be superimposed. 5308c XYZ1(3,NATM) represents the coordinates of MOL1 5309c XYZ2(3,NATM) represents the coordinates of MOL2 5310C A(3*3) represents the output matrix. 5311c T(3*3) represents the output vector. 5312c 5313C COMMON/SHIFS/ SCALR,SCALT 5314c DATA SCALR/60./,SCALT/1./ 5315C Note: There should be a scale of the shifts. i.e. Shift= scale * shifts 5316c Here SCALR is the rotation shiftscale 5317c Here SCALT is the translation shiftscale 5318c The initial and suggested SCALR is 60. 5319c The initial and suggested SCALt is 1. 5320C 5321C COMMON/RMS/ RMS,SMEAN,NREF,NREF1 5322C RMS is the final R.M.S between two molecules. 5323C SMEAN is the mean distance. 5324c NREF is the number of cycles of rotation refinement 5325c NREF1 is not useful 5326c They could be used to judge the SCALR and SCALS 5327c 5328C COMMON/IAT/ IAT1,IAT2,IAT3,IAT 5329C DATA SCALR/60./,IAT/0/ 5330C The routine uses the first three atoms to get the initial Eulerian 5331C angle. IAT1, IAT2 and IAT3 are the indexes of the three selected atoms. 5332c If IAT = 0, these atoms will be selected inside this routine. If not, 5333c the three atoms will be defined outside the routine through the common 5334c block. If the program does not work or you want special atoms, change 5335c IAT1,IAT2,IAT3 in COMMON/IAT/ and set IAT to 0 in order to select 5336c the three atoms in the main routine. 5337c 5338C NATM cannot be larger than NATM0 (currently 50000) or the parameter NATM0 5339C should be modified in this routine. 5340c 5341c The program use translation linear to least square refinement. Testing 5342c shows that it works quite well. The function of this routine is similar 5343c to SUPOS but the r.m.s. is less and the orientation might be more dangerous. 5344c 5345c 5346c by Guoguang Lu 5347C 27/01/1989 5348C 5349C 5350 COMMON/RMS/ RMS,SMEAN,NREF,NREF1 5351 COMMON/SHIFS/ SCALR,SCALS 5352 COMMON/IAT/ IAT1,IAT2,IAT3,IAT 5353c DATA SCALR/60./,SCALS/1./ 5354C DATA IAT/0/,IAT1/1/,IAT2/2/,IAT3/3/ 5355 PARAMETER (NATM0 = 50000) 5356c NATM0 is the largest number of atoms. 5357 DIMENSION XYZ1(3,NATM), XYZ2(3,NATM), A(3,3), T(3) 5358 DIMENSION X1(3,NATM0),X2(3,NATM0),VT(3*6*NATM0),DIS(3*NATM0) 5359 DIMENSION CEN1(3),CEN2(3) 5360 DIMENSION B1(3),B2(3),B3(3),T0(3) 5361c 5362C Number of atoms cannot be more than NATM0 or less than 3. 5363c 5364 IF (NATM.GT.NATM0) THEN 5365 WRITE(6,*) 'ERROR> Atom is more than ',NATM0 5366 STOP 5367 ELSE IF (NATM.LT.3) THEN 5368 STOP 'ERROR> Atom is less than 3.' 5369 END IF 5370 5371c Compute the initial matrix. 5372 CALL ORIEN(NATM,XYZ1,XYZ2,A) 5373C compute the gravity of mol1 and mol2 to CEN1 and CEN2 5374 CALL AVERG(3,NATM,XYZ1,CEN1) 5375 CALL AVERG(3,NATM,XYZ2,CEN2) 5376C T is the initial translation vector 5377 CALL ARRPS(3,1,CEN2,CEN1,T0) 5378c Change the origin to CEN1. 5379 CALL TMOVE(3,NATM,XYZ1,CEN1,-1.,X1) 5380 CALL TMOVE(3,NATM,XYZ2,CEN1,-1.,X2) 5381c 5382C Refine th1, th2, th3, tx, ty and tz 5383c 5384 CALL REFRT( NATM, X1, X2, A, T0, VT, DIS) 5385c 5386 WRITE(6,*) 5387 WRITE(6,*) 'R.M.S.' 5388 WRITE(6,*) ' natm' 5389 WRITE(6,*) 'SQRT( SIGMA (Distance(i)^2)/natm ) = ',RMS 5390 WRITE(6,*) ' i=1' 5391C 5392 WRITE(6,*) 5393 WRITE(6,*) 'Mean Distance:' 5394 WRITE(6,*) ' natm' 5395 WRITE(6,*) 'SIGMA (Distance(i)) /natm = ',SMEAN 5396 WRITE(6,*) ' i=1' 5397C 5398 WRITE(6,'(A)') '0' 5399 WRITE(6,'(A)') '0' 5400 WRITE(6,*) 'Mol1 is superposed to Mol2.' 5401 WRITE(6,*) 'The matrix and the vector are:' 5402 WRITE(6,*) 5403 WRITE(6,1) (A(1,J),J=1,3),CEN1(1),T0(1) 5404 WRITE(6,2) (A(2,J),J=1,3),CEN1(2),T0(2) 5405 WRITE(6,1) (A(3,J),J=1,3),CEN1(3),T0(3) 54061 FORMAT(1X,' (',3F10.6,' ) ( ',F8.3,' ) (' 5407 1 ,F8.3,' )') 54082 FORMAT(1X,' X2 = (',3F10.6,' ) * ( X1 -',F8.3,' ) + (' 5409 1 ,F8.3,' )') 5410C 5411 CALL MATMULT(3,3,3,1,A,CEN1,B3) 5412 CALL ARRPS(3,1,T0,B3,B3) 5413 CALL ARRAD(3,1,CEN1,B3,T) 5414C 5415 WRITE(6,'(A)') '0' 5416 WRITE(6,'(A)') '0' 5417 WRITE(6,*) 5418 WRITE(6,4) (A(1,J),J=1,3),T(1) 5419 WRITE(6,5) (A(2,J),J=1,3),T(2) 5420 WRITE(6,4) (A(3,J),J=1,3),T(3) 54214 FORMAT(1X,' (',3F10.6,' ) ( ) (',F8.3,' )') 54225 FORMAT(1X,' X2 = (',3F10.6,' ) * ( X1 ) + (',F8.3,' )') 5423C 5424 END 5425c 5426c 5427 SUBROUTINE SUPRIMP ( NATM, XYZ1, XYZ2, A, T ) 5428C DIMENSION XYZ1(3,NATM), XYZ2(3,NATM), A(3,3), T(3) 5429C A subroutine to give the superimpose matrix and vector from MOL1 to 5430c MOL2 i.e. XYZ2(3*NATM) = A(3*3) * XYZ1(3*NATM) + T(3) 5431C Where NATM represents the number of the atoms in each molecule 5432c which will be superimposed. 5433c XYZ1(3,NATM) represents the coordinates of MOL1 5434c XYZ2(3,NATM) represents the coordinates of MOL2 5435C A(3*3) represents the output matrix. 5436c T(3*3) represents the output vector. 5437c 5438C NATM cannot be larger than NATM0 (currently 50000) or the parameter NATM0 5439C should be modified in this routine. 5440c 5441C COMMON/SHIFS/ SCALR,SCALT 5442c DATA SCALR/1./,SCALT/1./ 5443C Note: There should be a scale of the shifts. i.e. Shift= scale * shifts 5444c Here SCALR is the rotation shiftscale 5445c Here SCALT is the translation shiftscale 5446c The initial and suggested SCALR is 1. 5447c The initial and suggested SCALt is 1. 5448C 5449C COMMON/RMS/ RMS,SMEAN,NREF,NREF1 5450C RMS is the final R.M.S between two molecules. 5451C SMEAN is the mean distance. 5452c NREF is the number of cycles of rotation refinement 5453c NREF1 is not used by this routine. 5454c They could be used to judge the SCALR and SCALS 5455c 5456C COMMON/IAT/ IAT1,IAT2,IAT3,IAT 5457C DATA SCALR/1./,IAT/0/ 5458C The routine use the first three atoms to get the initial Eulerian 5459C angle. IAT1, IAT2 and IAT3 are the indexes of the three select atoms. 5460c If IAT = 0, these atoms will be selected inside this routine. If not, 5461c the three atoms will be defined outside the routine through the common 5462c block. If the program does not work or you want special atoms, change 5463c IAT1,IAT2,IAT3 in COMMON/IAT/ and set IAT to 0 in order to select 5464c the three atoms in the main routine. Subroutine ORIEN performs this 5465c computation. 5466c 5467c Then the program use the subroutine REFORN to refine the orientation 5468c by vector method. The refinement equation is same with subroutine SUPOS. 5469C 5470c The program use translation linear and Eulerian non-linear least square 5471c refinement to refine both th1,th2,th3 and tx,ty,tz. Testing shows 5472c that it can give very low r.m.s., much less than any other method 5473c in most cases, however it is a little expensive in computer time. 5474c 5475c by Guoguang Lu 5476C 27/01/1989 5477C 5478 COMMON/RMS/ RMS,SMEAN,NREF,NREF1 5479 COMMON/SHIFS/ SCALR,SCALS 5480 COMMON/IAT/ IAT1,IAT2,IAT3,IAT 5481c DATA SCALR/1./,SCALS/1./,IAT/0/,IAT1/1/,IAT2/2/,IAT3/3/ 5482 PARAMETER (NATM0 = 50000) 5483c NATM0 is the largest number of atoms. 5484 DIMENSION XYZ1(3,NATM), XYZ2(3,NATM), A(3,3), T(3) 5485 DIMENSION X1(3,NATM0),X2(3,NATM0),VT(3*6*NATM0),DIS(3*NATM0) 5486 DIMENSION CEN1(3),CEN2(3) 5487 DIMENSION B1(3),B2(3),B3(3),T0(3) 5488c 5489c DEG=180./3.14159265359 5490c 5491C Number of atoms cannot be more than NATM0 or less than 3. 5492c 5493 IF (NATM.GT.50000) THEN 5494 WRITE(6,*) 'ERROR> Atom is more than ',NATM0 5495 STOP 5496 ELSE IF (NATM.LT.3) THEN 5497 STOP 'ERROR> Atom is less than 3.' 5498 END IF 5499 5500c Compute the initial matrix. 5501 CALL ORIEN(NATM,XYZ1,XYZ2,A) 5502c NREF1 = NREF 5503C 5504C Refine the orientation by vector method. 5505C 5506 CALL REFORN ( NATM, XYZ1, XYZ2, A, VT, DIS) 5507C 5508C compute the gravity of mol1 and mol2 to CEN1 and CEN2 5509 CALL AVERG(3,NATM,XYZ1,CEN1) 5510 CALL AVERG(3,NATM,XYZ2,CEN2) 5511C T is the initial translation vector 5512 CALL ARRPS(3,1,CEN2,CEN1,T0) 5513c Change the origin to CEN1. 5514 CALL TMOVE(3,NATM,XYZ1,CEN1,-1.,X1) 5515 CALL TMOVE(3,NATM,XYZ2,CEN1,-1.,X2) 5516c 5517C Refine th1, th2, th3, tx, ty and tz 5518c 5519 CALL REFRT( NATM, X1, X2, A, T0, VT, DIS) 5520c 5521 WRITE(6,*) 5522 WRITE(6,*) 'R.M.S.' 5523 WRITE(6,*) ' natm' 5524 WRITE(6,*) 'SQRT( SIGMA (Distance(i)^2)/natm ) = ',RMS 5525 WRITE(6,*) ' i=1' 5526C 5527 WRITE(6,*) 5528 WRITE(6,*) 'Mean Distance:' 5529 WRITE(6,*) ' natm' 5530 WRITE(6,*) 'SIGMA (Distance(i)) /natm = ',SMEAN 5531 WRITE(6,*) ' i=1' 5532C 5533 WRITE(6,*) 5534 WRITE(6,*) 'Mol1 is superposed to Mol2.' 5535 WRITE(6,*) 'The matrix and the vector are:' 5536 WRITE(6,*) 5537 WRITE(6,1) (A(1,J),J=1,3),CEN1(1),T0(1) 5538 WRITE(6,2) (A(2,J),J=1,3),CEN1(2),T0(2) 5539 WRITE(6,1) (A(3,J),J=1,3),CEN1(3),T0(3) 55401 FORMAT(1X,' (',3F10.6,' ) ( ',f10.5,' ) (' 5541 1 ,f10.5,' )') 55422 FORMAT(1X,' X2 = (',3F10.6,' ) * ( X1 -',f10.5,' ) + (' 5543 1 ,f10.5,' )') 5544C 5545 CALL MATMULT(3,3,3,1,A,CEN1,B3) 5546 CALL ARRPS(3,1,T0,B3,B3) 5547 CALL ARRAD(3,1,CEN1,B3,T) 5548C 5549 WRITE(6,*) 5550 WRITE(6,*) 5551 WRITE(6,4) (A(1,J),J=1,3),T(1) 5552 WRITE(6,5) (A(2,J),J=1,3),T(2) 5553 WRITE(6,4) (A(3,J),J=1,3),T(3) 55544 FORMAT(1X,' (',3F10.6,' ) ( ) (',f10.5,' )') 55555 FORMAT(1X,' X2 = (',3F10.6,' ) * ( X1 ) + (',f10.5,' )') 5556C 5557 END 5558c 5559c 5560 SUBROUTINE SUPRIMPFIN ( NATM, XYZ1, XYZ2, A, T ) 5561C DIMENSION XYZ1(3,NATM), XYZ2(3,NATM), A(3,3), T(3) 5562C A subroutine to give the superimpose matrix and vector from MOL1 to 5563c MOL2 i.e. XYZ2(3*NATM) = A(3*3) * XYZ1(3*NATM) + T(3) 5564C Where NATM represents the number of the atoms in each molecule 5565c which will be superimposed. 5566c XYZ1(3,NATM) represents the coordinates of MOL1 5567c XYZ2(3,NATM) represents the coordinates of MOL2 5568C A(3*3) represents the output matrix. 5569c T(3*3) represents the output vector. 5570c 5571C NATM cannot be larger than NATM0 (currently 50000) or the parameter NATM0 5572C should be modified in this routine. 5573c 5574C COMMON/SHIFS/ SCALR,SCALT 5575c DATA SCALR/1./,SCALT/1./ 5576C Note: There should be a scale of the shifts. i.e. Shift= scale * shifts 5577c Here SCALR is the rotation shiftscale 5578c Here SCALT is the translation shiftscale 5579c The initial and suggested SCALR is 1. 5580c The initial and suggested SCALt is 1. 5581C 5582C COMMON/RMS/ RMS,SMEAN,NREF,NREF1 5583C RMS is the final R.M.S between two molecules. 5584C SMEAN is the mean distance. 5585c NREF is the number of the cycles of the rotation refinement 5586c NREF1 is not used by this routine. 5587c They could be used to judge the SCALR,and SCALS 5588c 5589C COMMON/IAT/ IAT1,IAT2,IAT3,IAT 5590C DATA SCALR/1./,IAT/0/ 5591C The routine use the first three atoms to get the initial Eulerian 5592C angle. IAT1, IAT2 and IAT3 are the indexes of the three select atoms. 5593c If IAT = 0, these atoms will be selected inside this routine. If not, 5594c the three atoms will be defined outside the routine through the common 5595c block. If the program does not work or you want special atoms, change 5596c IAT1,IAT2,IAT3 in COMMON/IAT/ and give the IAT to 0 in order to select 5597c the three atoms in the main routine. Subroutine ORIEN is to proceed this 5598c computing 5599c 5600c Then the program uses the subroutine REFORN to refine the orientation 5601c by vector method. The refinement equation is same with subroutine SUPOS. 5602C 5603c The program use translation linear and Eulerian non-linear least square 5604c refinement to refine both th1,th2,th3 and tx,ty,tz. Testing shows 5605c that it can give very low r.m.s., much less than any other method 5606c in most cases, however it is a little expensive in computer time. 5607c 5608c by Guoguang Lu 5609C 27/01/1989 5610C 5611 COMMON/RMS/ RMS,SMEAN,NREF,NREF1 5612 COMMON/SHIFS/ SCALR,SCALS 5613 COMMON/IAT/ IAT1,IAT2,IAT3,IAT 5614c DATA SCALR/1./,SCALS/1./,IAT/0/,IAT1/1/,IAT2/2/,IAT3/3/ 5615 PARAMETER (NATM0 = 50000) 5616c NATM0 is the maximum number of atoms. 5617 DIMENSION XYZ1(3,NATM), XYZ2(3,NATM), A(3,3), T(3) 5618 DIMENSION X1(3,NATM0),X2(3,NATM0),VT(3*6*NATM0),DIS(3*NATM0) 5619 DIMENSION CEN1(3),CEN2(3) 5620 DIMENSION B1(3),B2(3),B3(3),T0(3) 5621c 5622c DEG=180./3.14159265359 5623c 5624C Number of atoms cannot be more than NATM0 or less than 3. 5625c 5626 IF (NATM.GT.50000) THEN 5627 WRITE(6,*) 'ERROR> Atom is more than ',NATM0 5628 STOP 5629 ELSE IF (NATM.LT.3) THEN 5630 STOP 'ERROR> Atom is less than 3.' 5631 END IF 5632 5633c Compute the initial matrix. 5634 CALL ORIEN(NATM,XYZ1,XYZ2,A) 5635c NREF1 = NREF 5636C 5637C Refine the orientation by vector method. 5638cc 5639 CALL REFORNFIN ( NATM, XYZ1, XYZ2, A, VT, DIS) 5640 NREFORN = NREF 5641C 5642C compute the gravity of mol1 and mol2 to CEN1 and CEN2 5643 CALL AVERG(3,NATM,XYZ1,CEN1) 5644 CALL AVERG(3,NATM,XYZ2,CEN2) 5645C T is the initial translation vector 5646 CALL ARRPS(3,1,CEN2,CEN1,T0) 5647c Change the origin to CEN1. 5648 CALL TMOVE(3,NATM,XYZ1,CEN1,-1.,X1) 5649 CALL TMOVE(3,NATM,XYZ2,CEN1,-1.,X2) 5650c 5651C Refine th1, th2, th3, tx, ty and tz 5652c 5653 CALL REFRTFIN ( NATM, X1, X2, A, T0, VT, DIS ) 5654 NREFRT = NREF1 5655 CALL REFRTFIN1 ( NATM, X1, X2, A, T0, VT, DIS, VT(1+2*6*NATM) ) 5656 NREFRT1 = NREF1 5657c 5658 WRITE(6,*) 'Final fit cycle>>>', NREFORN, NREFRT, NREFRT1 5659 WRITE(6,*) 5660 WRITE(6,*) 'R.M.S.' 5661 WRITE(6,*) ' natm' 5662 WRITE(6,*) 'SQRT( SIGMA (Distance(i)^2)/natm ) = ',RMS 5663 WRITE(6,*) ' i=1' 5664C 5665 WRITE(6,*) 5666 WRITE(6,*) 'Mean Distance:' 5667 WRITE(6,*) ' natm' 5668 WRITE(6,*) 'SIGMA (Distance(i)) /natm = ',SMEAN 5669 WRITE(6,*) ' i=1' 5670C 5671 WRITE(6,*) 5672 WRITE(6,*) 'Mol1 is superposed to Mol2.' 5673 WRITE(6,*) 'The matrix and the vector are:' 5674 WRITE(6,*) 5675 WRITE(6,1) (A(1,J),J=1,3),CEN1(1),T0(1) 5676 WRITE(6,2) (A(2,J),J=1,3),CEN1(2),T0(2) 5677 WRITE(6,1) (A(3,J),J=1,3),CEN1(3),T0(3) 56781 FORMAT(1X,' (',3F10.6,' ) ( ',f10.5,' ) (' 5679 1 ,f10.5,' )') 56802 FORMAT(1X,' X2 = (',3F10.6,' ) * ( X1 -',f10.5,' ) + (' 5681 1 ,f10.5,' )') 5682C 5683 CALL MATMULT(3,3,3,1,A,CEN1,B3) 5684 CALL ARRPS(3,1,T0,B3,B3) 5685 CALL ARRAD(3,1,CEN1,B3,T) 5686C 5687 WRITE(6,*) 5688 WRITE(6,*) 5689 WRITE(6,4) (A(1,J),J=1,3),T(1) 5690 WRITE(6,5) (A(2,J),J=1,3),T(2) 5691 WRITE(6,4) (A(3,J),J=1,3),T(3) 56924 FORMAT(1X,' (',3F10.6,' ) ( ) (',f10.5,' )') 56935 FORMAT(1X,' X2 = (',3F10.6,' ) * ( X1 ) + (',f10.5,' )') 5694C 5695 END 5696c 5697c 5698 subroutine lgg_symlib(nunit,file,nsp,nsym,nrot,sym,latnam) 5699c In this subroutine the user gives a space group number, the subroutine reads the 5700c ccp4 symmetry library, then it outputs the symmetry operations 5701c Input 5702c file--- symmetry library file name 5703c nsp --- space group number 5704c Output 5705c nsymm --- Total number of symmetry operations in this spacegroup 5706c nrot ---- only rotation symmetric operations in this spacegroup 5707c symm ---- a 3*4*nsym symmetric operation c matrix. 5708c latnam --- *14 character spacegroup name 5709c in common symdump, if dump is true, dump all the information 5710c maxop is the maximum number of symmetry operations 5711c libary format 5712c155 6 2 R32 5713c X,Y,Z * -Y,X-Y,Z * Y-X,-X,Z 5714c Y,X,-Z * -X,Y-X,-Z * X-Y,-Y,-Z 5715c X+1/3,Y+2/3,Z+2/3 * -Y+1/3,X-Y+2/3,Z+2/3 * Y-X+1/3,-X+2/3,Z+2/3 5716c Y+1/3,X+2/3,-Z+2/3 * -X+1/3,Y-X+2/3,-Z+2/3 * X-Y+1/3,-Y+2/3,-Z+2/3 5717c X+2/3,Y+1/3,Z+1/3 * -Y+2/3,X-Y+1/3,Z+1/3 * Y-X+2/3,-X+1/3,Z+1/3 5718c Y+2/3,X+1/3,-Z+1/3 * -X+2/3,Y-X+1/3,-Z+1/3 * X-Y+2/3,-Y+1/3,-Z+1/3 5719 parameter (maxop=96) 5720 character file*80,latnam*14 5721 real sym(3,4,maxop) 5722 integer*4 nunit,nsp,nsym,nrot 5723 logical*1 dump 5724c COMMON/SYMDUMP/ DUMP 5725c 5726 character*80 key,str*40 5727c 5728c... data statements. Separate declaration and init required for f2c 5729 data dump /.false./ 5730c 5731c WRITE(6,*) 'nsp',nsp 5732 nsym=0 5733 if (nunit.eq.0) nunit=25 5734 if (file(1:1).eq.' ') call spstrunct(file) 5735 if (file(1:1).eq.' ') file='SYMOP' 5736 call ccpdpn(nunit,file(1:lenstr(file)),'READONLY','F',0,0) 573710 read(nunit,'(a)',end=220) key 5738c il = lenstr(key) 5739c i1 = index(key(1:il),'PG') 5740c if (i1.eq.0) goto 10 5741c call redstrin(12,i1-1,key,npar) 5742 read(key(1:12),*,err=10,end=10) isp,iall,irot 5743c WRITE(6,*) key(1:lenstr(key)),isp,nsp,latnam 5744 if (isp.ne.nsp) goto 10 5745c i1 = index(key(1:il),'PG') 5746 read(key(12:80),'(a)') latnam 5747c call spstrunct(latnam) 5748 if (latnam(1:1).eq.' ') latnam(1:13) = latnam(2:14) 5749 if (isp.ne.nsp) goto 10 5750 write(6,*) 'Space Group >>> ',Latnam,isp 5751 do i = 1, iall 5752 read(nunit,'(a)') key 5753 il = lenstr(key) 5754 key(il+1:il+1)='*' 5755 iend = 0 5756 if (dump) write(6,'(1x,a)') key(1:il) 575720 ist = iend+1 5758 iend = index(key(ist:il+1),'*')+IST-1 5759 str = ' ' 5760 str(1:iend-ist) = key(ist:iend-1) 5761 NSYM=NSYM+1 5762 IF (MATSYM(SYM(1,1,nsym),STR,ICOL).EQ.0) GOTO 502 5763 WRITE(6,*) 'Error in symop after colunm ',ICOL 5764 write(6,*) key 5765 write(6,*) str 5766 stop 'Check you file SYMOP' 5767502 continue 5768 if (dump) then 5769 write(6,'(4f8.4)') ((sym(j1,j2,nsym),j2=1,4),j1=1,3) 5770 end if 5771 if (iend.lt.il+1) goto 20 5772 if (i.eq.irot) nrot = nsym 5773 end do 5774 write(6,45) nsym,nrot 577545 format(1x,'Symmetric operation ----',6x, 5776 1 'Total: ',i3,6x,'Rotation:',i3) 5777 close (unit=nunit) 5778 return 5779210 WRITE(6,*) 'Error opening the SYMOP file ',FILE(1:LENSTR(FILE)) 5780220 WRITE(6,*) 'Space group',nsp,' was not found in SYMOP file' 5781 nsym = 0 5782 nrot = 0 5783 latnam = ' ' 5784 STOP 'Check your SYMOP file' 5785 END 5786c 5787c 5788 SUBROUTINE TMOVE(M,NATM,XIN,T,CONST,XOUT) 5789C XOUT = XIN + CONST * T 5790C where XIN is input M*NATM-dimensional array. 5791C XOUT N is output M*NATM-dimensional array. 5792c T is a M-dimensional vector. 5793c CONST is a constant. 5794C DIMENSION XIN(M,NATM),XOUT(M,NATM),T(M),B(100) 5795c 5796 DIMENSION XIN(M,NATM),XOUT(M,NATM),T(M),B(100) 5797 CALL ARRMC(M,1,T,CONST,B) 5798 DO I = 1, NATM 5799 CALL ARRAD(M,1,XIN(1,I),B,XOUT(1,I)) 5800 END DO 5801 END 5802c 5803c 5804c------------------------------------------------------------ 5805 subroutine up(txt,len) 5806c 5807c convert character string to upper case 5808c 5809 character*(*) txt 5810 character*80 save 5811 character*26 ualpha,lalpha 5812 data ualpha /'ABCDEFGHIJKLMNOPQRSTUVWXYZ'/ 5813 data lalpha /'abcdefghijklmnopqrstuvwxyz'/ 5814 5815 do 9000 i=1,len 5816 if((txt(i:i).ge.'a').and.(txt(i:i).le.'z')) then 5817 match = index(lalpha,txt(i:i)) 5818 save(i:i) = ualpha(match:match) 5819 else 5820 save(i:i) = txt(i:i) 5821 end if 5822c end do 58239000 continue 5824 5825 txt = save 5826 return 5827 end 5828c 5829c 5830 SUBROUTINE VECCRSMLT(A1,A2,A) 5831 DIMENSION A1(3),A2(3),A(3) 5832 A(1)=VLDIM2(A1(2),A1(3),A2(2),A2(3)) 5833 A(2)=VLDIM2(A1(3),A1(1),A2(3),A2(1)) 5834 A(3)=VLDIM2(A1(1),A1(2),A2(1),A2(2)) 5835 END 5836c 5837c 5838 FUNCTION VEM(N,V) 5839 DIMENSION V(N) 5840 C=0 5841 DO 10 I=1,N 584210 C=C+V(I)*V(I) 5843 VEM=SQRT(C) 5844 RETURN 5845 END 5846c 5847c 5848 FUNCTION VLDIM3(AT) 5849C A function to calculate the modulus of a 3*3-dimension matrix. 5850c VLDIM3 = | AT3*3 | 5851 dimension at(3,3) 5852 dimension B1(3) 5853C 5854 call veccrsmlt(at(1,1),at(1,2),b1) 5855 vldim3 = poimult(3,3,b1,at(1,3)) 5856 end 5857c 5858c 5859 FUNCTION VLDIM2(A11,A12,A21,A22) 5860 VLDIM2=A11*A22-A21*A12 5861 END 5862c 5863c 5864 REAL FUNCTION COSD(ANGINDEG) 5865 REAL ANGINDEG 5866 COSD = COS(ANGINDEG*3.1415926/180.) 5867 END 5868c 5869c 5870 REAL FUNCTION SIND(ANGINDEG) 5871 REAL ANGINDEG 5872 SIND = SIN(ANGINDEG*3.1415926/180.) 5873 END 5874c 5875c 5876 REAL FUNCTION TAND(ANGINDEG) 5877 REAL ANGINDEG 5878 TAND = TAN(ANGINDEG*3.1415926/180.) 5879 END 5880c 5881c 5882 REAL FUNCTION ACOSD(VAL) 5883 REAL VAL 5884 ACOSD = ACOS(VAL)*180./3.1415926 5885 END 5886c 5887c 5888 REAL FUNCTION ASIND(VAL) 5889 REAL VAL 5890 ASIND = ASIN(VAL)*180./3.1415926 5891 END 5892c 5893c 5894 REAL FUNCTION ATAND(VAL) 5895 REAL VAL 5896 ATAND = ATAN(VAL)*180./3.1415926 5897 END 5898