C C This code is distributed under the terms and conditions of the C CCP4 licence agreement as `Part ii)' software. See the conditions C in the CCP4 manual for a copyright statement. C C********************** C Program ABS ** C********************** C C Purpose: to determine the absolute configuration of the heavy atom C substructure or polar space group using anomalous scattering C PROGRAM ABSC PARAMETER(MAXREF=90000) PARAMETER(MAXDATA=4194304) COMPLEX FF(MAXREF),DATA(MAXDATA) INTEGER*2 IH(MAXREF),IK(MAXREF),IL(MAXREF) LOGICAL EOF, LVALMS CHARACTER*4 KEY, CCVAL(20),KEYWORD(4) CHARACTER*400 LINE CHARACTER*80 TITLE INTEGER NTOK, IBEG(20), IEND(20), ITYP(20), IDEC(20) REAL FVAL(20) DIMENSION ADATA1(7),X(480,64),Y(480,64),Z(480,64),OCU(480) INTEGER IPOINT1(7) CHARACTER*30 LAB1(7),LSUSRJ(7) CHARACTER*1 CTYP1(7) REAL RSYM(4,4,64),HKL(3) DATA LAB1/'H','K','L','F','SIGF','DANO','SIGDANO'/ DATA KEYWORD/'TITL','ATOM','LABI','RESO'/ DATA CTYP1/'H','H','H','F','Q','D','Q'/ DATA IPOINT1/3*-1,4*0/ DATA OCU/480*1.0/ CALL CCPRCS (6, 'ABS', '$Date: 2002/11/04 16:41:21 $') CALL CCPFYP C Set default - currently NAN CALL QNAN (VAL_MAGIC) WRITE(6,*) WRITE(6,*)' Program: ABS ' WRITE(6,*)' Last update: 10/2/2001' WRITE(6,*)' Authors: Quan HAO' WRITE(6,*) TWOPI=2.*3.1415926 JJ = 0 1 LINE = ' ' NTOK = 20 CALL PARSER (KEY,LINE,IBEG,IEND,ITYP,FVAL,CCVAL, + IDEC,NTOK,EOF,.FALSE.) IF (EOF) GOTO 10 CALL CCPUPC(KEY) DO I = 1 , 4 IF (KEY.EQ.KEYWORD(I)) THEN GOTO (3,4,5,6) I END IF END DO CALL CCPERR(1,'----------- illegal keyword ------------') 3 TITLE = LINE(1:80) WRITE (6,*) TITLE WRITE (6,*) GOTO 1 4 JJ = JJ + 1 X(JJ,1) = FVAL(2) Y(JJ,1) = FVAL(3) Z(JJ,1) = FVAL(4) IF(NTOK.GE.5) OCU(JJ) = FVAL(5) write(6,800)jj, x(jj,1),y(jj,1),z(jj,1),ocu(jj) GOTO 1 5 CALL MTZINI ITOK = 2 CALL LKYSET(LAB1,7,LSUSRJ,IPOINT1,ITOK,NTOK,LINE, + IBEG,IEND) CALL LKYIN(1,LAB1,7,NTOK,LINE,IBEG,IEND) CALL LROPEN (1,'HKLIN',1,IFAIL) IF (IFAIL.EQ.1) + CALL CCPERR(1,'----------- FILE failure ------------') CALL LRSYMM(1,NSYM,RSYM) CALL LRASSN (1,LAB1,7,IPOINT1,CTYP1) GOTO 1 6 resolim = fval(2) goto 1 10 DO 14 J = 1, JJ DO 12 N = 2, NSYM X(J,N) = RSYM(1,1,N)*X(J,1) + RSYM(1,2,N)*Y(J,1) + 1 RSYM(1,3,N)*Z(J,1) + RSYM(1,4,N) Y(J,N) = RSYM(2,1,N)*X(J,1) + RSYM(2,2,N)*Y(J,1) + 1 RSYM(2,3,N)*Z(J,1) + RSYM(2,4,N) Z(J,N) = RSYM(3,1,N)*X(J,1) + RSYM(3,2,N)*Y(J,1) + 1 RSYM(3,3,N)*Z(J,1) + RSYM(3,4,N) 12 CONTINUE 14 CONTINUE NREF = 0 MAXH = 0 MAXK = 0 MAXL = 0 cen = 0. 20 CALL LRREFF(1,RESOL,ADATA1,EOF) IF (EOF) GOTO 50 CALL IS_MAGIC (VAL_MAGIC,adata1(4),LVALMS) IF (lvalms) goto 20 if (adata1(4).lt.0.) GOTO 20 CALL IS_MAGIC (VAL_MAGIC,adata1(6),LVALMS) IF (lvalms) adata1(6) = 0.0 if(1./sqrt(resol).lt.resolim) goto 20 NREF = NREF + 1 SUM_R = 0.0 SUM_I = 0.0 DO 24 J = 1, JJ DO 22 N = 1, NSYM ARG=TWOPI*(ADATA1(1)*X(J,N)+ . ADATA1(2)*Y(J,N)+ . ADATA1(3)*Z(J,N)) SUM_R=SUM_R+OCU(J)*COS(ARG) SUM_I=SUM_I+OCU(J)*SIN(ARG) 22 CONTINUE 24 CONTINUE DO 38 N = 1, NSYM DO 36 I = 1, 3 HKL(I) = RSYM(1,I,N)*ADATA1(1)+RSYM(2,I,N)*ADATA1(2) 1 +RSYM(3,I,N)*ADATA1(3) 36 CONTINUE MAXH = MAX(NINT(ABS(HKL(1))),MAXH) MAXK = MAX(NINT(ABS(HKL(2))),MAXK) MAXL = MAX(NINT(ABS(HKL(3))),MAXL) 38 CONTINUE IH(NREF) = NINT(ADATA1(1)) IK(NREF) = NINT(ADATA1(2)) IL(NREF) = NINT(ADATA1(3)) C + 90 DEGREES FF(NREF) = CMPLX(-SUM_I*ADATA1(4)*ADATA1(6), 1 SUM_R*ADATA1(4)*ADATA1(6)) c --- phase of ha substructure phi_ha=atan2(sum_i,sum_r) c write(6,*)ih(nref),ik(nref),il(nref),sum_r,sum_i cen = cen + abs(sin(phi_ha)) GOTO 20 50 DO I = 8, 0, -1 NTEMP = 2**I IF(MAXH.GE.NTEMP) GOTO 60 END DO 60 MAXH = NTEMP*4 DO I = 8, 0, -1 NTEMP = 2**I IF(MAXK.GE.NTEMP) GOTO 70 END DO 70 MAXK = NTEMP*4 DO I = 8, 0, -1 NTEMP = 2**I IF(MAXL.GE.NTEMP) GOTO 80 END DO 80 MAXL = NTEMP*4 cen = cen/float(nref) WRITE(6,1000)MAXH,MAXK,MAXL IF(MAXH*MAXK*MAXL.GT.MAXDATA) + CALL CCPERR(1,'*out of allocated memory, use RESO to limit + the number of reflections**') write (6,900) cen CALL FFTN(NSYM,RSYM,MAXH,MAXK,MAXL,NREF,IH,IK,IL,FF,DATA) CALL CCPERR(0,'Normal termination') 800 FORMAT('Heavy atom', i4,' : X =', 1 F7.4,', Y =', F7.4, ', Z =', F7.4, 2 '; occ =', F7.4,/) 900 format('<|sin(phase_ha_substructure)|> =', f8.5) 1000 FORMAT('Sampling intervals on xyz : ',3I8,/) END C********************************************************************** SUBROUTINE FFTN(NSYM,RSYM,MAXH,MAXK,MAXL,NREF,IH,IK,IL,FF,DATA) PARAMETER(MAXREF=90000) COMPLEX FF(MAXREF), DATA(0:MAXH-1,0:MAXK-1,0:MAXL-1) INTEGER*2 IH(MAXREF),IK(MAXREF),IL(MAXREF) INTEGER NN(3) REAL RSYM(4,4,64) REAL D2R PARAMETER (D2R = 3.1415927 / 180. ) DO I = 0, MAXH-1 DO J = 0, MAXK-1 DO K = 0, MAXL-1 DATA(I,J,K) = 0. END DO END DO END DO DO 100 I = 1, NREF DO 80 N = 1, NSYM DO 70 ISIGN = 1, -1, -2 PHIS = 360.*ISIGN*(RSYM(1,4,N)*IH(I) 1 +RSYM(2,4,N)*IK(I)+RSYM(3,4,N)*IL(I)) CPHIS = COS(PHIS*D2R) SPHIS = SIN(PHIS*D2R) JH = ISIGN*NINT(RSYM(1,1,N)*IH(I)+RSYM(2,1,N)*IK(I) 1 +RSYM(3,1,N)*IL(I)) JH = MOD(JH+MAXH, MAXH) JK = ISIGN*NINT(RSYM(1,2,N)*IH(I)+RSYM(2,2,N)*IK(I) 1 +RSYM(3,2,N)*IL(I)) JK = MOD(JK+MAXK, MAXK) JL = ISIGN*NINT(RSYM(1,3,N)*IH(I)+RSYM(2,3,N)*IK(I) 1 +RSYM(3,3,N)*IL(I)) JL = MOD(JL+MAXL, MAXL) A = REAL(FF(I)) B = ISIGN*AIMAG(FF(I)) DATA(JH,JK,JL) = CMPLX(A*CPHIS-B*SPHIS, A*SPHIS+B*CPHIS) 70 CONTINUE 80 CONTINUE 100 CONTINUE NN(1) = MAXH NN(2) = MAXK NN(3) = MAXL CALL FOURN(DATA,NN,3,-1) RHO3 = 0.0 RHOABS3 = 0.0 DO 180 I = 0, MAXH-1 DO 170 J = 0, MAXK-1 DO 160 K = 0, MAXL-1 RHO3 = RHO3 + REAL(DATA(I,J,K))**3 RHOABS3 = RHOABS3 + CABS(DATA(I,J,K))**3 160 CONTINUE 170 CONTINUE 180 CONTINUE C = RHO3/RHOABS3 IF(C.GT.0.) THEN WRITE(6,1010) C ELSE WRITE(6,1020) C END IF RETURN 1010 FORMAT(' C = ', F10.4,' , *correct configuration*') 1020 FORMAT(' C = ', F10.4,' , *incorrect configuration*') END C*********************************************************************** SUBROUTINE fourn(data,nn,ndim,isign) INTEGER isign,ndim,nn(ndim) REAL data(*) INTEGER i1,i2,i2rev,i3,i3rev,ibit,idim,ifp1,ifp2,ip1,ip2,ip3,k1, *k2,n,nprev,nrem,ntot REAL tempi,tempr DOUBLE PRECISION theta,wi,wpi,wpr,wr,wtemp ntot=1 do 11 idim=1,ndim ntot=ntot*nn(idim) 11 continue nprev=1 do 18 idim=1,ndim n=nn(idim) nrem=ntot/(n*nprev) ip1=2*nprev ip2=ip1*n ip3=ip2*nrem i2rev=1 do 14 i2=1,ip2,ip1 if(i2.lt.i2rev)then do 13 i1=i2,i2+ip1-2,2 do 12 i3=i1,ip3,ip2 i3rev=i2rev+i3-i2 tempr=data(i3) tempi=data(i3+1) data(i3)=data(i3rev) data(i3+1)=data(i3rev+1) data(i3rev)=tempr data(i3rev+1)=tempi 12 continue 13 continue endif ibit=ip2/2 1 if ((ibit.ge.ip1).and.(i2rev.gt.ibit)) then i2rev=i2rev-ibit ibit=ibit/2 goto 1 endif i2rev=i2rev+ibit 14 continue ifp1=ip1 2 if(ifp1.lt.ip2)then ifp2=2*ifp1 theta=isign*6.28318530717959d0/(ifp2/ip1) wpr=-2.d0*sin(0.5d0*theta)**2 wpi=sin(theta) wr=1.d0 wi=0.d0 do 17 i3=1,ifp1,ip1 do 16 i1=i3,i3+ip1-2,2 do 15 i2=i1,ip3,ifp2 k1=i2 k2=k1+ifp1 tempr=sngl(wr)*data(k2)-sngl(wi)*data(k2+1) tempi=sngl(wr)*data(k2+1)+sngl(wi)*data(k2) data(k2)=data(k1)-tempr data(k2+1)=data(k1+1)-tempi data(k1)=data(k1)+tempr data(k1+1)=data(k1+1)+tempi 15 continue 16 continue wtemp=wr wr=wr*wpr-wi*wpi+wr wi=wi*wpr+wtemp*wpi+wi 17 continue ifp1=ifp2 goto 2 endif nprev=n*nprev 18 continue return END