CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SAPI C A direct-method program to locate heavy atom positions. C Modified by Q. Hao (1,2) C Original Authors: Y.X.Gu, C.D.Zheng, J.X.Yao, J.Z.Qian & H.F.Fan (2) C Email: qh22@cornell.edu or fan@aphy.iphy.ac.cn C (1) Cornell University, USA C (2) Institute of Physics, Chinese Academy of Sciences, Beijing,China CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PROGRAM SAPI CHARACTER*132 AA CHARACTER*8 CDATE,CTIME C GET INITIAL TIME CALL CCPDAT(CDATE) CALL UTIME(CTIME) WRITE(6,910) CDATE,CTIME 910 FORMAT(53X,A8,3X,A8) CALL CCPFYP c call ccprcs (6, 'SAPI', '$Date: 2001/02/12 09:25:49 $') WRITE (6,100) CALL CCPDPN(10,'KEYWORD.TM','SCRATCH','F',80,0) c OPEN(10,FILE='KEYWORD.TM',FORM='FORMATTED',STATUS='UNKNOWN') 1 READ (5,10,END=2)AA WRITE (10,10)AA GOTO 1 2 CONTINUE CALL PREPAR CALL PHASE CALL EXFFT CALL SEARCH4 10 FORMAT(A132) 100 FORMAT(///,' ****** SAPI ******',/ & ' Version 1.1, 10/10/2001',// & 'A direct-method program to locate heavy atom positions.'// & 'Authors: YX Gu, JX Yao, CD Zheng, HF Fan (1) & Q Hao(2)',/ & '(1) Institute of Physics, Chinese Academy of Sciences',/ & '(2) Cornell University, USA.',/ & 'Email: qh22@cornell.edu or fan@aphy.iphy.ac.cn'//) CALL CCPERR(0,'Normal termination') END ************************************************************************ * * * PPPPPPPP RRRRRRRR EEEEEEEEE PPPPPPPP A RRRRRRRR * * PPPPPPPPP RRRRRRRRR EEEEEEEEE PPPPPPPPP AAA RRRRRRRRR * * PP PP RR RR EE PP PP AA AA RR RR * * PP PP RR RR EE PP PP AA AA RR RR * * PPPPPPPPP RRRRRRRRR EEEEEEE PPPPPPPPP AA AA RRRRRRRRR * * PPPPPPPP RRRRRRRR EEEEEEE PPPPPPPP AA AA RRRRRRRR * * PP RR RR EE PP AAAAAAAAA RR RR * * PP RR RR EE PP AAAAAAAAA RR RR * * PP RR RR EEEEEEEEE PP AA AA RR RR * * PP RR RR EEEEEEEEE PP AA AA RR RR * * * * PROGRAM FOR PRELIMINARY PROCESSING OF THE INPUT DATA * * ** AN EXTENSIVE MODIFICATION OF THE PROGRAM 'NORMAL' OF MULTAN-80 ** * * PC VERSION 1998 * ************************************************************************ SUBROUTINE PREPAR C INCLUDE 'FLIB.FI' COMMON/ATMFACTOR/AL(8),AS(8),BL(8),BS(8),CL(8),CS(8),DL(8),DS(8), 1 EL(8),NW(8),NO(8),NK,NAT,F(9) COMMON/ATMGROUP/NINF(10),NAG(10),X(5000),Y(5000),Z(5000),NZ(5000), 1 Q(5000),U11(5000),U22(5000),U33(5000),U23(5000), 2 U13(5000),U12(5000) COMMON/REFLXIN/LH(600),LK(600),LL(600),FO(600),ID(600),RHO(600), 1 SIG(600),XKP(600) COMMON/REFLXOUT/IX(30000),EX(30000),FX(30000),SCMK(30000) COMMON/RESCALING/ISC,IBGR,MG,MIG(8),SCL,BT(9),SC(9),IP(8,3,5), 1 SCAL(8),RSTT COMMON/SCATFACTOR/GIS(142),GIW(142),NGP,NDIFF,ISOL,X0(3) COMMON/SINETABLE/SINT(450) COMMON/SYMMETRY/IS(3,3,24),TS(3,24),NSYM,PTS,KSYS,ICENT,LATT, 1 IAVE,EXTI COMMON/STATISTICS/VST(10,5),NST(5),ZT(25,5),EE(10),MULT,IND,NZR, 1 TMUL C UNITS FOR INPUT/OUTPUT, TITLE, FLAGS COMMON/UNIT1/ ITLE(80),LIST,PI,KCURV,IPAT,IAPA,MOV,MFP,NAU COMMON/WILSON/FLGW(30,9),FLGD(30,9),AVR(30,9),DCV(50,9),SLOPE(9), 1 FLGK(9),DEL(9),KS(9),EW(600),ED(600),EDP(600) COMMON/XXX/P(6),CX(9),NREF,NB,RHOMAX,MM,EMAX,EN,MZ,ER,ISIM,RHOCUT, 1 RHOMIN,RHOLOW,ENG(8),ERG(8),KNOWN,IPATH,VOLUME CHARACTER NGS(26),KX(10),ITERM(4),NTLE(80),ITLE,KSP,KP,KM,KEQ,KC DATA KX/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/ DATA KSP/1H /,KP/1H+/,KM/1H-/,KEQ/1H=/,KSC/1H;/ DATA ITERM/1HH,1HK,1HL,1HN/ C SET UP INITIAL VALUES, READ PROGRAM PARAMETERS PI=4.0*ATAN(1.0) NREF=0 RHOMAX=0.0 RHOMIN=1.0 DTOR=PI/180.0 C SET UP SIN/COS TABLE DO 50 I=1,450 SINT(I)=SIN(DTOR*FLOAT(I-1)) 50 CONTINUE CALL CCPDPN(2,'FCOEF.TM','SCRATCH','U',80,0) CALL CCPDPN(8,'SCRA8.TM','SCRATCH','U',80,0) c OPEN(2,FILE='FCOEF.TM',FORM='UNFORMATTED',STATUS='UNKNOWN') C OPEN(5,FILE='KEYWORD.TM',FORM='FORMATTED',STATUS='UNKNOWN') C OPEN(6,FILE='PREPARE.OUT',FORM='FORMATTED',STATUS='UNKNOWN') c OPEN(8,FILE='SCRA8.TM',FORM='UNFORMATTED',STATUS='UNKNOWN') REWIND 10 C READ TITLE READ (10,100) NTLE 100 FORMAT(80A1) IL=80 IR=80 DO 103 I=1,80 IF (NTLE(I).NE.KSP.AND.IL.EQ.80) IL=I-1 IF (NTLE(80-I+1).NE.KSP.AND.IR.EQ.80) IR=I-1 103 CONTINUE IM=(IL+IR)/2 DO 105 I=1,80 IF (I.LE.IM) ITLE(I)=KSP IF (I.GT.IM) ITLE(I)=NTLE(IL+I-IM) 105 CONTINUE c-qh IF (IERR.EQ.1) GO TO 900 CALL INPUT_P(JUMP,NON,MISS,KMARK,NTOTAL,BTT) IF(NON.NE.1) GO TO 300 STOP' --- PREPARE COMPLETED ---' C READ REFLECTION DATA FROM DATA FILE 300 CALL DATAIN(MH,MK,ML) ANAT=FLOAT(NAT)/(PTS*FLOAT((ICENT+1)*NSYM)) WRITE (6,370) ANAT 370 FORMAT(/1X,36HNUMBER OF ATOMS IN ASYMMETRIC UNIT =,F8.2) NASU=INT(ANAT+0.5) C CHECK IF WEIGHTED OR DIFFERENCE FOURIER IS REQUIRED NDIFF=0 DO 400 I=1,NGP IF (NINF(I).EQ.6.OR.NINF(I).EQ.7) GO TO 410 400 CONTINUE IF (MM.GT.0) WRITE (6,405) MM,MZ 405 FORMAT(1X,21HOUTPUT FOR PHASE - ,I5,13H LARGEST E'S,3X, 1 'AND',I5,14H SMALLEST E'S) GO TO 420 410 NDIFF=1 REWIND 8 WRITE (2) NDIFF WRITE (2) ITLE,ICENT,LATT,NSYM,((TS(I,J),(IS(K,I,J),K=1,3), 1 I=1,3),J=1,NSYM),(NDIFF,I=1,3),(CX(I),I=1,6),NASU,RHOCUT 2 ,RHOLOW EN=0.0 ER=0.0 DO 415 I=1,8 ENG(I)=0.0 ERG(I)=0.0 415 CONTINUE SC(1)=SQRT(SCAL(1)) IF (SC(1).LT.0.01) SC(1)=1.0 420 MKR=0 IF (NDIFF.EQ.1) CALL RFAC(ICENT,NINF) IF (NDIFF.EQ.1) GO TO 900 C CALCULATE NUMBER OF REFLEXIONS TO PASS TO PHASE IF (MM.GT.0) GO TO 422 MN = INT(16.0*ANAT+150.5) IF (ICENT.EQ.1) MN = MN + 50 IF (NSYM.EQ.1) MN = MN + 50 MM = MIN0(MN,15000) IF (MN.LE.500) MM=MN+MIN0(500-MN,100)/2 MM=MAX0(MM,200) C 422 MM=MIN0(MM,NREF-100) MM=MIN0(MM,NREF-100) 422 IF (NREF.LE.500.AND.EN.EQ.1.2) EN=0.9 IF (NREF.LE.300.AND.MZ.EQ.100) MZ=50 IF (ISC.NE.2) GO TO 635 C SPECIAL RESCALING MKR=MKR+1 432 WRITE (6,500) MG 500 FORMAT(/1X,13X,41HREFLECTIONS ARE RESCALED ACCORDING TO & ,I3,13H INDEX GROUPS) IF (MKR.EQ.1) WRITE (6,510) 510 FORMAT(28X,23H(SPECIFIED BY THE USER)) C PRINT OUT THE INDEX RELATIONSHIPS DO 591 MI=1,MG MJ=0 DO 530 J=1,3 DO 520 K=1,3 IF (IP(MI,J,K).EQ.0) GO TO 520 MJ=MJ+1 GO TO 530 520 CONTINUE 530 CONTINUE IF (MJ.EQ.0) GO TO 591 DO 580 J=1,MJ DO 554 II=1,26 NGS(II)=KSP 554 CONTINUE I=0 II=2 DO 570 K=1,5 IF (IP(MI,J,K).NE.0) GO TO 555 IF (K.NE.5) GO TO 560 IF (IP(MI,J,4).NE.0) GO TO 560 555 KK=IABS(IP(MI,J,K)) IF (I.EQ.0) GO TO 556 IF (IP(MI,J,K).GT.0) NGS(II)=KP 556 IF (IP(MI,J,K).LT.0) NGS(II)=KM I=1 IF (NGS(II).NE.KSP) II=II+2 IF (KK.EQ.1.AND.K.NE.5) GO TO 558 IF (KK.LT.10) GOTO 557 NGS(II)=KX(INT(KK/10)+1) II=II+1 KK=MOD(KK,10) 557 NGS(II)=KX(KK+1) II=II+2 558 IF (K.EQ.5) GO TO 560 NGS(II)=ITERM(K) II=II+2 560 IF (K.NE.3) GO TO 570 NGS(II)=KEQ I=0 II=II+2 570 CONTINUE 573 FORMAT(/25X,5HGROUP,I3,5X,26A1) 574 FORMAT(38X,26A1) IF (J.EQ.1) WRITE (6,573) MI,NGS IF (J.NE.1) WRITE (6,574) NGS 580 CONTINUE 591 CONTINUE 635 IF (NB.EQ.0) NB=8.0*ALOG10(0.05*FLOAT(MAX0(NREF,100))+0.5) C MAXIMUM OF 30 POINTS ON WILSON PLOT IF (NB.GT.30) NB=30 RHOCUT=AMIN1(RHOCUT,RHOMAX) RHOLOW=AMAX1(RHOLOW,RHOMIN) ANGMAX=0.5*SQRT(1.0/RHOMAX) ANGLHI=0.5*SQRT(1.0/RHOCUT) ANGLOW=0.5*SQRT(1.0/RHOLOW) WRITE (6,640) NREF,RHOMAX,ANGMAX,ANGLOW,ANGLHI,NB 640 FORMAT(/1X,'TOTAL NUMBER OF REFLEXION (IN DATA-FILE) =',I6/ 1 1X,37HMAXIMUM (SIN(THETA)/LAMBDA)**2 =,F9.4,' (',F5.2, 2 ') ANGSTROM'/1X,'LIMIT OF RESOLUTION FOR PHASE/EXFFT =', 3 F7.2,' to',F7.2,' ANGSTROM'/ 4 1X,33HNUMBER OF POINTS ON WILSON PLOT =,I3) SC(1)=1.0 IF (JUMP.GE.0) GO TO 650 C OBTAIN SUMS FOR WILSON PLOT AND FIT LEAST SQUARES STRAIGHT LINE CALL WILSONN(MH,MK,ML,MISS,BTT) IF (ISC.EQ.2.OR.NAU.EQ.1) GO TO 645 CALL AUTOGP(LATT,EN) IF (ISC.EQ.2) GO TO 432 C PLOT WILSON AND DEBYE CURVES AND LEAST SQUARES STRAIGHT LINE 645 CALL GRAPH(NGP) 650 DO 655 I=1,8 BT(I)=2.0*BT(I) IF (BT(I).LT.-199.9) BT(I)=BT(1) IF (SCAL(I).LT.0.00001) SCAL(I)=SCAL(1) 655 CONTINUE IF (IPAT.NE.0) GO TO 750 C CALCULATE SCALE FACTORS FOR APPROPRIATE REFLEXION GROUPS CALL RESCA(KSYS,JUMP) IF (NINF(1).EQ.4) NGP=0 C OUTPUT FOR PHASE CALL CCPDPN(1,'PHASDT.TM','SCRATCH','F',80,0) c OPEN(1,FILE='PHASDT.TM',FORM='FORMATTED',STATUS='UNKNOWN') WRITE (1,660) ITLE,(CX(I),I=1,6),NW,NO,ICENT,LATT,NSYM,NGP,NASU 660 FORMAT(80A1/6F10.5/16I5/5I5) WRITE (1,670) ((TS(I,J),(IS(K,I,J),K=1,3),I=1,3),J=1,NSYM) 670 FORMAT(3(F11.8,3I3)) IF (NGP.EQ.0) GO TO 700 NF=0 DO 690 I=1,NGP NS=NF+1 NF=NF+NAG(I) WRITE (1,680) NINF(I),NAG(I),(NZ(J),X(J),Y(J),Z(J),J=NS,NF) 680 FORMAT(2I5/(I5,3F10.6)) 690 CONTINUE C CALCULATE FINAL E'S AND OUTPUT REFLEXION STATISTICS 700 CALL ECAL IF (NGP.LE.0) GO TO 720 NF=0 DO 710 IGP=1,NGP NS=NF+1 NF=NF+NAG(IGP) IF (NINF(IGP).LT.5) GO TO 710 IF (NINF(IGP).EQ.5) CALL RECYC(NS,NF,KMARK) GO TO 720 710 CONTINUE 720 CONTINUE c CLOSE (1) GO TO 900 750 NDIFF=2 REWIND 2 WRITE (2) NDIFF WRITE (2) ITLE,ICENT,LATT,NSYM,((TS(I,J),(IS(K,I,J),K=1,3), 1 I=1,3),J=1,NSYM),(NDIFF,I=1,3),(CX(I),I=1,6),NASU,RHOCUT 2 ,RHOLOW CALL PATT 900 CONTINUE c CLOSE (2) close (8) c CLOSE (8,STATUS='DELETE') C OUTPUT FOR SEARCH REWIND 11 c CALL CCPDPN(11,'SEARCH.TM','UNKNOWN','F',80,0) CALL CCPDPN(38,'SEARKW.TM','SCRATCH','F',80,0) c OPEN(UNIT=11,FILE='SEARCH.TM',FORM='FORMATTED',STATUS='UNKNOWN') c OPEN(UNIT= 8,FILE='SEARKW.TM',FORM='FORMATTED',STATUS='UNKNOWN') DO 950 I=1,NTOTAL READ (11,100) NTLE WRITE (38,100) NTLE 950 CONTINUE CLOSE (11) c CLOSE (38) END C----------------------------------------------------------------------- SUBROUTINE ATMCOEF C CALCULATE SPHERICAL SCATTERING FACTOR COEFFICIENTS C ATOMIC SCATTERING FACTOR COEFFICIENTS FOR 98 ATOM TYPES C F = AL * EXP(-AS * RHO) + BL * EXP(-BS * RHO) C + CL * EXP(-CS * RHO) + DL * EXP(-DS * RHO) + EL C COMMON/ATMFACTOR/AL(8),AS(8),BL(8),BS(8),CL(8),CS(8),DL(8),DS(8), 1 EL(8),NW(8),NO(8),NK,NAT,F(9) DIMENSION ALT(98),AST(98),BLT(98),BST(98),ELT(98),N1(98) DIMENSION CLT(98),CST(98),DLT(98),DST(98) DATA N1/ 8,805,1209,205, 2, 3, 14, 15, 6,1405,1401,1307,112,1909, 1 16, 19, 312,118, 11, 301,1903,2009, 22, 318,1314,605,315,1409, 2 321,2614, 701,705,119,1905,218,1118,1802,1918,25,2618,1402,1315, 3 2003,1821,1808,1604, 107,304,914,1914,1902,2005, 9,2405,319,201, 4 1201,305,1618,1404,1613,1913,521,704,2002,425,815,518,2013,2502, 5 1221,806,2001, 23,1805,1519,918,1620,121,807,2012,1602,209,1615, 6 120,1814,618,1801, 103,2008,1601, 21,1416,1621,113,313, 211,306/ DATA ALT/0.489918,0.873400, 1.128200, 1.591900, 2.054500, 9 2.310000, 12.212600, 3.048500, 3.539200, 3.955300, 9 4.762600, 5.420400, 6.420200, 6.291500, 6.434500, 9 6.905300, 11.460400, 7.484500, 8.218600, 8.626600, 9 9.189000, 9.759500, 10.297100, 10.640600, 11.281900, 9 11.769500, 12.284100, 12.837600, 13.338000, 14.074300, 9 15.235400, 16.081600, 16.672300, 17.000601, 17.178900, 9 17.355499, 17.178400, 17.566299, 17.775999, 17.876499, 9 17.614201, 3.702500, 19.130100, 19.267401, 19.295700, 9 19.331900, 19.280800, 19.221399, 19.162399, 19.188900, 9 19.641800, 19.964399, 20.147200, 20.293301, 20.389200, 9 20.336100, 20.577999, 21.167101, 22.044001, 22.684500, 9 23.340500, 24.004200, 24.627399, 25.070900, 25.897600, 9 26.507000, 26.904900, 27.656300, 28.181900, 28.664101, 9 28.947599, 29.143999, 29.202400, 29.081800, 28.762100, 9 28.189400, 27.304899, 27.005899, 16.881901, 20.680901, 9 27.544600, 31.061701, 33.368900, 34.672600, 35.316299, 9 35.563099, 35.929901, 35.763000, 35.659698, 35.564499, 9 35.884701, 36.022800, 36.187401, 36.525398, 36.670601, 9 36.648800, 36.788101, 36.918499/ DATA AST/20.659300,9.103700,3.954600, 43.642700, 23.218500, 9 20.843901, 0.005700, 13.277100, 10.282500, 8.404200, 9 3.285000, 2.827500, 3.038700, 2.438600, 1.906700, 9 1.467900, 0.010400, 0.907200, 12.794900, 10.442100, 9 9.021300, 7.850800, 6.865700, 6.103800, 5.340900, 9 4.761100, 4.279100, 3.878500, 3.582800, 3.265500, 9 3.066900, 2.850900, 2.634500, 2.409800, 2.172300, 9 1.938400, 1.788800, 1.556400, 1.402900, 1.276180, 9 1.188650, 0.277200, 0.864132, 0.808520, 0.751536, 9 0.698655, 0.644600, 0.594600, 0.547600, 5.830300, 9 5.303400, 4.817420, 4.347000, 3.928200, 3.569000, 9 3.216000, 2.948170, 2.812190, 2.773930, 2.662480, 9 2.562700, 2.472740, 2.387900, 2.253410, 2.242560, 9 2.180200, 2.070510, 2.073560, 2.028590, 1.988900, 9 1.901820, 1.832620, 1.773330, 1.720290, 1.671910, 9 1.629030, 1.592790, 1.512930, 0.461100, 0.545000, 9 0.655150, 0.690200, 0.704000, 0.700999, 0.685870, 9 0.663100, 0.646453, 0.616341, 0.589092, 0.563359, 9 0.547751, 0.529300, 0.511929, 0.499384, 0.483629, 9 0.465154, 0.451018, 0.437533/ DATA BLT/0.262003,0.630900, 0.750800, 1.127800, 1.332600, 9 1.020000, 3.132200, 2.286800, 2.641200, 3.112500, 9 3.173600, 2.173500, 1.900200, 3.035300, 4.179100, 9 5.203400, 7.196400, 6.772300, 7.439800, 7.387300, 9 7.367900, 7.355800, 7.351100, 7.353700, 7.357300, 9 7.357300, 7.340900, 7.292000, 7.167600, 7.031800, 9 6.700600, 6.374700, 6.070100, 5.819600, 5.235800, 9 6.728600, 9.643500, 9.818400, 10.294600, 10.948000, 9 12.014400, 17.235600, 11.094800, 12.918200, 14.350100, 9 15.501700, 16.688499, 17.644400, 18.559601, 19.100500, 9 19.045500, 19.013800, 18.994900, 19.029800, 19.106199, 9 19.297001, 19.599001, 19.769501, 19.669701, 19.684700, 9 19.609501, 19.425800, 19.088600, 19.079800, 18.218500, 9 17.638300, 17.294001, 16.428499, 15.885100, 15.434500, 9 15.220800, 15.172600, 15.229300, 15.430000, 15.718900, 9 16.155001, 16.729601, 17.763901, 18.591299, 19.041700, 9 19.158400, 13.063700, 12.951000, 15.473300, 19.021099, 9 21.281601, 23.054701, 22.906401, 23.103201, 23.421900, 9 23.294800, 23.412800, 23.596399, 23.808300, 24.099199, 9 24.409599, 24.773600, 25.199499/ DATA BST/7.740390,3.356800, 1.052400, 1.862300, 1.021000, 9 10.207500, 9.893300, 5.701100, 4.294400, 3.426200, 9 8.842200, 79.261101, 0.742600, 32.333698, 27.157000, 9 22.215099, 1.166200, 14.840700, 0.774800, 0.659900, 9 0.572900, 0.500000, 0.438500, 0.392000, 0.343200, 9 0.307200, 0.278400, 0.256500, 0.247000, 0.233300, 9 0.241200, 0.251600, 0.264700, 0.272600, 16.579599, 9 16.562300, 17.315100, 14.098800, 12.800600, 11.916000, 9 11.766000, 1.095800, 8.144870, 8.434670, 8.217580, 9 7.989290, 7.472600, 6.908900, 6.377600, 0.503100, 9 0.460700, 0.420885, 0.381400, 0.344000, 0.310700, 9 0.275600, 0.244475, 0.226836, 0.222087, 0.210628, 9 0.202088, 0.196451, 0.194200, 0.181951, 0.196143, 9 0.202172, 0.197940, 0.223545, 0.238849, 0.257119, 9 9.985190, 9.599900, 9.370460, 9.225900, 9.092270, 9 8.979480, 8.865530, 8.811740, 8.621600, 8.448400, 9 8.707510, 2.357600, 2.923800, 3.550780, 3.974580, 9 4.069100, 4.176190, 3.871350, 3.651550, 3.462040, 9 3.415190, 3.325300, 3.253960, 3.263710, 3.206470, 9 3.089970, 3.046190, 3.007750/ DATA CLT/0.196767,0.311200, 0.617500, 0.539100, 1.097900, 9 1.588600, 2.012500, 1.546300, 1.517000, 1.454600, 9 1.267400, 1.226900, 1.593600, 1.989100, 1.780000, 9 1.437900, 6.255600, 0.653900, 1.051900, 1.589900, 9 1.640900, 1.699100, 2.070300, 3.324000, 3.019300, 9 3.522200, 4.003400, 4.443800, 5.615800, 5.165200, 9 4.359100, 3.706800, 3.431300, 3.973100, 5.637700, 9 5.549300, 5.139900, 5.422000, 5.726290, 5.417320, 9 4.041830, 12.887600, 4.649010, 4.863370, 4.734250, 9 5.295370, 4.804500, 4.461000, 4.294800, 4.458500, 9 5.037100, 6.144870, 7.513800, 8.976700, 10.662000, 9 10.888000, 11.372700, 11.851300, 12.385600, 12.774000, 9 13.123500, 13.439600, 13.760300, 13.851800, 14.316700, 9 14.559600, 14.558300, 14.977900, 15.154200, 15.308700, 9 15.100000, 14.758600, 14.513500, 14.432700, 14.556400, 9 14.930500, 15.611500, 15.713100, 25.558201, 21.657499, 9 15.538000, 18.441999, 16.587700, 13.113800, 9.498870, 9 8.003700, 12.143900, 12.473900, 12.597700, 12.747300, 9 14.189100, 14.949100, 15.640200, 16.770700, 17.341499, 9 17.399000, 17.891899, 18.331699/ DATA CST/49.551899,22.927601,85.390503,103.483002, 60.349800, 9 0.568700, 28.997499, 0.323900, 0.261500, 0.230600, 9 0.313600, 0.380800, 31.547199, 0.678500, 0.526000, 9 0.253600, 18.519400, 43.898300, 213.186996, 85.748398, 9 136.108002, 35.633801, 26.893801, 20.262600, 17.867399, 9 15.353500, 13.535900, 12.176300, 11.396600, 10.316300, 9 10.780500, 11.446800, 12.947900, 15.237200, 0.260900, 9 0.226100, 0.274800, 0.166400, 0.125599, 0.117622, 9 0.204785, 11.004000, 21.570700, 24.799700, 25.874901, 9 25.205200, 24.660500, 24.700800, 25.849899, 26.890900, 9 27.907400, 28.528400, 27.766001, 26.465900, 24.387899, 9 20.207300, 18.772600, 17.608299, 16.766899, 15.885000, 9 15.100900, 14.399600, 13.754600, 12.933100, 12.664800, 9 12.189900, 11.440700, 11.360400, 10.997500, 10.664700, 9 0.261033, 0.275116, 0.295977, 0.321703, 0.350500, 9 0.382661, 0.417916, 0.424593, 1.482600, 1.572900, 9 1.963470, 8.618000, 8.793700, 9.556420, 11.382400, 9 14.042200, 23.105200, 19.988701, 18.599001, 17.830900, 9 16.923500, 16.092699, 15.362200, 14.945500, 14.313600, 9 13.434600, 12.894600, 12.404400/ DATA DLT/0.049879,0.178000, 0.465300, 0.702900, 0.706800, 9 0.865000, 1.166300, 0.867000, 1.024300, 1.125100, 9 1.112800, 2.307300, 1.964600, 1.541000, 1.490800, 9 1.586300, 1.645500, 1.644200, 0.865900, 1.021100, 9 1.468000, 1.902100, 2.057100, 1.492200, 2.244100, 9 2.304500, 2.348800, 2.380000, 1.673500, 2.410000, 9 2.962300, 3.683000, 4.277900, 4.354300, 3.985100, 9 3.537500, 1.529200, 2.669400, 3.265880, 3.657210, 9 3.533460, 3.742900, 2.712630, 1.567560, 1.289180, 9 0.605844, 1.046300, 1.602900, 2.039600, 2.466300, 9 2.682700, 2.523900, 2.273500, 1.990000, 1.495300, 9 2.695900, 3.287190, 3.330490, 2.824280, 2.851370, 9 2.875160, 2.896040, 2.922700, 3.545450, 2.953540, 9 2.965770, 3.638370, 2.982330, 2.987060, 2.989630, 9 3.716010, 4.300130, 4.764920, 5.119820, 5.441740, 9 5.675890, 5.833770, 5.783700, 5.860000, 5.967600, 9 5.525930, 5.969600, 6.469200, 7.025880, 7.425180, 9 7.443300, 2.112530, 3.210970, 4.086550, 4.807030, 9 4.172870, 4.188000, 4.185500, 3.479470, 3.493310, 9 4.216650, 4.232840, 4.243910/ DATA DST/2.201590,0.982100,168.261002, 0.542000, 0.140300, 9 51.651199, 0.582600, 32.908901, 26.147600, 21.718399, 9 129.423996, 7.193700, 85.088600, 81.693703, 68.164497, 9 56.172001, 47.778400, 33.392899, 41.684101, 178.436996, 9 51.353100, 116.105003, 102.477997, 98.739899, 83.754303, 9 76.880501, 71.169197, 66.342102, 64.812599, 58.709702, 9 61.413502, 54.762501, 47.797199, 43.816299, 41.432800, 9 39.397202, 164.934006, 132.376007, 104.353996, 87.662697, 9 69.795700, 61.658401, 86.847198, 94.292801, 98.606201, 9 76.898598, 99.815598, 87.482498, 92.802902, 83.957100, 9 75.282501, 70.840302, 66.877602, 64.265800, 213.904007, 9 167.201996, 133.123993, 127.112999, 143.643997, 137.903000, 9 132.720993, 128.007004, 123.174004, 101.398003, 115.362000, 9 111.874001, 92.656601, 105.703003, 102.960999, 100.417000, 9 84.329803, 72.028999, 63.364399, 57.056000, 52.086102, 9 48.164700, 45.001099, 38.610298, 36.395599, 38.324600, 9 45.814899, 47.257900, 48.009300, 47.004501, 45.471500, 9 44.247299, 150.645004, 142.324997, 117.019997, 99.172203, 9 105.250999, 100.612999, 97.490799, 105.980003, 102.273003, 9 88.483398, 86.002998, 83.788101/ DATA ELT/0.001305,0.006400, 0.037700, 0.038500, -0.193200, 9 0.215600, -11.529000, 0.250800, 0.277600, 0.351500, 9 0.676000, 0.858400, 1.115100, 1.140700, 1.114900, 9 0.866900, -9.557400, 1.444500, 1.422800, 1.375100, 9 1.332900, 1.280700, 1.219900, 1.183200, 1.089600, 9 1.036900, 1.011800, 1.034100, 1.191000, 1.304100, 9 1.718900, 2.131300, 2.531000, 2.840900, 2.955700, 9 2.825000, 3.487300, 2.506400, 1.912130, 2.069290, 9 3.755910, 4.387500, 5.404280, 5.378740, 5.328000, 9 5.265930, 5.179000, 5.069400, 4.939100, 4.782100, 9 4.590900, 4.352000, 4.071200, 3.711800, 3.335200, 9 2.773100, 2.146780, 1.862640, 2.058300, 1.984860, 9 2.028760, 2.209630, 2.574500, 2.419600, 3.583240, 9 4.297280, 4.567960, 5.920460, 6.756210, 7.566720, 9 7.976280, 8.581540, 9.243540, 9.887500, 10.472000, 9 11.000500, 11.472200, 11.688300, 12.065800, 12.608900, 9 13.174600, 13.411800, 13.578200, 13.677000, 13.710800, 9 13.690500, 13.724700, 13.621100, 13.526600, 13.431400, 9 13.428700, 13.396600, 13.357300, 13.381200, 13.359200, 9 13.288700, 13.275400, 13.267400/ DO 150 I=1,NK C CHECK ATOM TYPE DO 120 J=1,98 IF (NW(I).NE.N1(J)) GO TO 120 NO(I)=J TEST=AL(I)+BL(I)+CL(I)+ABS(DL(I))+ABS(EL(I)) C TEST IF THE PARAMETERS HAVE BEEN INPUT BY USER IN A KEYWORD FILE IF (TEST.GT.0.001) GO TO 120 AS(I)=AST(J) AL(I)=ALT(J) BS(I)=BST(J) BL(I)=BLT(J) CL(I)=CLT(J) CS(I)=CST(J) DL(I)=DLT(J) DS(I)=DST(J) EL(I)=ELT(J) 120 CONTINUE 150 CONTINUE RETURN END ************************************************************************ * * * A U U TTTTTTT OOOOO GGGGG PPPPPP * * A A U U T O O G G P P * * A A U U T O O G P P * * A A U U T O O G GGGG PPPPPP * * AAAAAAA U U T O O G G P * * A A U U T O O G G P * * A A UUUUU T OOOOO GGGGGG P * * * * SEARCH FOR THE PSEUDO-SYSTEMATIC EXTINCTION RULE AND * * PREPARE FOR AUTOMATIC GROUPING OF THE REFLECTIONS * * VERSION 1998 * ************************************************************************ SUBROUTINE AUTOGP(LATT,EMIN) COMMON/REFLXIN/LH(600),LK(600),LL(600),FO(600),ID(600),RHO(600), 1 SIG(600),XKP(600) COMMON/REFLXOUT/INX(3,30000),EX(30000) COMMON/RESCALING/ISC,IBGR,MG,MIG(8),SCL,BT(9),SC(9),IP(8,3,5), 1 SCAL(8) COMMON/WILSON/FLGW(30,9),FLGD(30,9),AVR(30,9),DCV(50,9),SLOPE(9), 1 FLGK(9),DEL(9),KS(9),EW(600),ED(600),EDP(600) DIMENSION EDX(3,256),EDX1(3,256),MDX(3),ICX(3),IN(3), 1 IP6(6,5),ILSYM(3) WRITE (6,10) 10 FORMAT(//9X,45H*** AUTOMATIC SEARCH FOR PSEUDO - SYSTEMATIC , & 14HEXTINCTION ***) DO 40 I=1,3 IN(I)=0 ICX(I)=1 ILSYM(I)=2 DO 30 J=1,256 EDX(I,J)=0.0 EDX1(I,J)=0.0 30 CONTINUE 40 CONTINUE DO 50 I=1,6 DO 50 J=1,5 IP6(I,J)=0 50 CONTINUE C STORE THE SYSTEMATIC EXTINCTION RULES DUE TO LATTICE CENTERING IF (LATT.EQ.1) GO TO 200 K=LATT-1 GO TO (100,110,120,200,200,150), K 100 ILSYM(1)=1 GO TO 200 110 ILSYM(2)=1 GO TO 200 120 ILSYM(3)=1 GO TO 200 150 DO 160 I=1,3 ILSYM(I)=3 160 CONTINUE C READ REFLECTION FILE, STORE SEPARATIVELY THE ABSOLUTE VALUES OF C INDICES H, K, L AND THE CUMULATED E**2 VALUES 200 M=MG+1 NB=0 REWIND 8 210 READ (8) LH,LK,LL,FO,ID,EW,ED,RHO DO 220 I=1,600 IF (FO(I).LT.0.0) GO TO 222 ED(I)=SC(1)*ED(I)*EXP(2.0*BT(1)*RHO(I)) IF (SQRT(ED(I)).LT.EMIN) GO TO 220 C PICK UP REFLECTIONS WITH E GREATER THAN EMIN IH=IABS(LH(I)) IK=IABS(LK(I)) IL=IABS(LL(I)) IF (IH.NE.0) EDX(1,IH)=EDX(1,IH)+ED(I) IF (IK.NE.0) EDX(2,IK)=EDX(2,IK)+ED(I) IF (IL.NE.0) EDX(3,IL)=EDX(3,IL)+ED(I) C STORE UP TO 15000 STRONG REFLECTIONS FOR LATER USE NB=NB+1 INX(1,NB)=LH(I) INX(2,NB)=LK(I) INX(3,NB)=LL(I) EX(NB)=ED(I) IF (NB.EQ.15000) GO TO 222 220 CONTINUE GO TO 210 222 WRITE (6,224) NB 224 FORMAT(24X,I6,2X,'LARGEST E**2 WERE USED'/) DO 230 I=1,3 DO 228 J=1,256 EDX1(I,J)=EDX(I,J) 228 CONTINUE 230 CONTINUE C FIND THE ONE-INDEX-RELATIONS FOR THE 'STRONG' RELECTION GROUP NRLN=0 CALL GCD(EDX,ICX,INDEX,SCL,IP,IP6,NRLN) IF (INDEX.EQ.1) GO TO 250 DO 240 I=1,3 IF (ICX(I).EQ.1) GO TO 240 NRLN=NRLN+1 IP6(NRLN,I)=1 IP6(NRLN,4)=ICX(I) ICX(I)=1 240 CONTINUE C FIND THE THREE INDEX-DIFFERENCES FOR EACH PAIR OF REFLECTIONS, C STORE ONE INDEX-DIFFERENCE WHEN THE OTHER TWO EQUAL ZERO 250 DO 265 I=1,3 DO 260 J=1,256 EDX(I,J)=0.0 260 CONTINUE 265 CONTINUE DO 450 I=1,NB N=I+1 DO 400 J=N,NB DO 310 K=1,3 MDX(K)=IABS(INX(K,I)-INX(K,J)) 310 CONTINUE DO 320 K=1,3 IF (MDX(K).EQ.0) GO TO 320 KSW=0 DO 314 KK=1,3 IF (KK.EQ.K) GO TO 314 IF (MDX(KK).NE.0) GO TO 314 KSW=KSW+1 IF (KSW.NE.2) GO TO 314 MX=MDX(K) EDX(K,MX)=EDX(K,MX)+EX(I)+EX(J) 314 CONTINUE 320 CONTINUE 400 CONTINUE 450 CONTINUE WRITE (6,456) 456 FORMAT(/1X,' H SIGMA K SIGMA L SIGMA',3X, & 'DELTA H SIGMA DELTA K SIGMA DELTA L SIGMA'/ & 1X,3(5X,4HE**2),1X,3(11X,4HE**2)/) WRITE (6,458)(I,EDX1(1,I),I,EDX1(2,I),I,EDX1(3,I),I,EDX(1,I), & I,EDX(2,I),I,EDX(3,I),I=1,16) 458 FORMAT((1X,3(I3,F6.1),2X,3(5X,I3,F7.1))) CALL GCD(EDX,ICX,INDEX,SCL,IP,IP6,NRLN) C REJECT INDEX-DIFFERENCE RELATIONS DUE TO LATTICE CENTERING IF (LATT.EQ.1) GO TO 510 DO 500 I=1,3 IF (ICX(I).EQ.1.OR.ICX(I).NE.ILSYM(I)) GO TO 500 ICX(I)=1 INDEX=INDEX-1 500 CONTINUE 510 IF (INDEX.GT.2) GO TO 680 IF (NRLN.GT.0) GO TO 1105 GO TO 1200 C FIND THE RELATIONS INVOLVING TWO INDICES 680 NRS=NRLN NRS2=NRLN DO 760 I=1,2 N=I+1 DO 740 J=N,3 DO 700 J1=1,3 IN(J1)=0 700 CONTINUE CALL COMMUL(ICX(I),ICX(J),1,ICOM) DO 720 K=2,3 IS=(-1)**K IN(I)=ICOM/ICX(I) IN(J)=IS*ICOM/ICX(J) CALL RCHECK(IN,ICOM,NRLN,JUMP,NB,IP6) IF (NRLN-NRS.EQ.2) GO TO 770 IF (JUMP.EQ.1) GO TO 740 IF (ICOM.EQ.2) GO TO 740 720 CONTINUE 740 CONTINUE 760 CONTINUE 770 NRS2=NRLN IF (INDEX.LT.4) GO TO 1000 C RELATION INVOLVING THREE INDICES CALL COMMUL(ICX(1),ICX(2),ICX(3),ICOM) IN(3)=ICOM/ICX(3) DO 800 I=2,3 IS=(-1)**I IN(1)=IS*ICOM/ICX(1) DO 790 K=2,3 KKS=(-1)**K IN(2)=KKS*ICOM/ICX(2) CALL RCHECK(IN,ICOM,NRLN,JUMP,NB,IP6) IF (JUMP.EQ.1) GO TO 1000 790 CONTINUE 800 CONTINUE 1000 IF (NRLN.EQ.NRS) GO TO 1105 C CHECK FOR THE INDEPENDENCY OF THE INDEX RELATIONS NRLNT=NRLN NRS1=NRS+1 IF (NRLN.EQ.NRS2.OR.NRS2.EQ.NRS) GO TO 1045 C SIMULTANEOUS EXISTENCE OF 2- AND 3-INDEX RELATIONS IS NOT ALLOWED DO 1020 I=NRS1,NRS2 IF (IP6(I,5).GT.IP6(NRLN,5)) GO TO 1020 C DELETE 3-INDEX RELATION DO 1010 J=1,5 IP6(NRLN,J)=0 1010 CONTINUE NRLN=NRLN-1 GO TO 1045 1020 CONTINUE C DELETE 2-INDEX RELATIONS DO 1040 I=NRS1,NRS2 DO 1030 J=1,5 IP6(I,J)=0 1030 CONTINUE NRLN=NRLN-1 1040 CONTINUE 1045 IF (NRS.EQ.0) GO TO 1105 C DELETE 1-INDEX RELATION IF IT IS INCLUDED IN ANOTHER RELATION OR C IN THE COMBINATION OF OTHER RELATIONS NR0=NRS DO 1100 I=1,NRS J=NRS1-I DO 1050 K=1,3 IF (IP6(J,K).EQ.1) GO TO 1052 1050 CONTINUE 1052 DO 1090 L=NRS1,NRLNT IF (IP6(L,4).EQ.0.OR.IP6(L,K).EQ.0) GO TO 1090 II=1 DO 1054 M=1,4 IF (M.EQ.K) GO TO 1054 IN(II)=IP6(L,M) II=II+1 1054 CONTINUE IC=IP6(L,K) CALL GCD3(IN,IC,JUMP) IF (JUMP.EQ.1) GO TO 1066 DO 1064 II=1,4 IP6(J,II)=0 1064 CONTINUE NR0=NR0-1 1066 IF (NR0.LT.2) GO TO 1090 DO 1070 II=1,3 IN(II)=0 IF (NRLNT.GT.NRS2.AND.IP6(NRLNT,4).NE.0) IN(II)=1 1070 CONTINUE II=1 DO 1080 M=1,3 IF (M.EQ.J) GO TO 1080 DO 1078 MM=1,NRS IF (IP6(MM,4).EQ.0) GO TO 1078 IF (MM.EQ.J) GO TO 1078 IF (IP6(MM,M).EQ.0.OR.IP6(L,M).EQ.0) GO TO 1078 IN(II)=IP6(MM,4)*IP6(L,MM) II=II+1 1078 CONTINUE 1080 CONTINUE IF (II.EQ.1) GO TO 1090 IN(II)=IP6(L,4) IC=IP6(L,K) CALL GCD3(IN,IC,JUMP) IF (JUMP.EQ.1) GO TO 1100 DO 1088 II=1,4 IP6(J,II)=0 1088 CONTINUE NR0=NR0-1 1090 CONTINUE 1100 CONTINUE 1105 NRLN=0 DO 1110 I=1,6 IF (IP6(I,4).EQ.0) GO TO 1110 NRLN=NRLN+1 DO 1120 J=1,4 IP(1,NRLN,J)=IP6(I,J) IP6(I,J)=0 1120 CONTINUE 1110 CONTINUE IF (NRLN.GT.0) GO TO 1500 1200 WRITE (6,1210) 1210 FORMAT(/1X,12X,45H*** NO PSEUDO-TRANSLATIONAL SYMMETRY HAS BEEN, & 10H FOUND ***/) RETURN C SET UP INDEX RELATIONS FOR THE WEAK GROUPS 1500 MG=1 DO 1510 I=1,NRLN MG=MG*IP(1,I,4) 1510 CONTINUE IF (MG.LE.8) GO TO 1520 MG=2 GO TO 1620 1520 DO 1600 I=2,MG DO 1580 J=1,NRLN DO 1560 K=1,4 IP(I,J,K)=IP(1,J,K) 1560 CONTINUE IF (NRLN-J.NE.2) GO TO 1565 IP(I,J,5)=MOD((I-1)/(IP(1,NRLN,4)*IP(1,NRLN-1,4)),IP(1,J,4)) GO TO 1580 1565 IF (NRLN-J.EQ.1) IP(I,J,5)=MOD((I-1)/IP(1,NRLN,4),IP(1,J,4)) IF (NRLN-J.EQ.0) IP(I,J,5)=MOD(I-1,IP(1,J,4)) 1580 CONTINUE 1600 CONTINUE 1620 ISC=2 IBGR=1 C CHECK FOR NULL 'WEAK' GROUPS II=0 MG1=MG DO 1700 NI=2,MG I=NI-II ESIG=0 REWIND 8 1630 READ (8) LH,LK,LL,FO,ID,EW,ED DO 1660 N=1,600 IF (FO(N).LT.0.0) GO TO 1665 DO 1650 J=1,3 IF (IP(I,J,4).EQ.0) GO TO 1650 IF (MOD(IP(I,J,1)*LH(N)+IP(I,J,2)*LK(N)+IP(I,J,3)*LL(N)-IP(I,J,5) & ,IP(I,J,4)).NE.0) GO TO 1660 1650 CONTINUE ESIG=ESIG+ED(N) 1660 CONTINUE N=N-1 1665 IF (ESIG.GT.0.000001) GO TO 1700 IF (FO(N).GE.0.0) GO TO 1630 MG1=MG1-1 DO 1690 K=I,MG1 DO 1680 J=1,3 DO 1670 JJ=1,5 IP(K,J,JJ)=IP(K+1,J,JJ) 1670 CONTINUE 1680 CONTINUE 1690 CONTINUE II=1+II 1700 CONTINUE MG=MG1 CALL SCRWRT WRITE (6,1800) 1800 FORMAT(//6X,48H*** PSEUDO-TRANSLATIONAL SYMMETRY HAS BEEN FOUND, & 19H BY THE PROGRAM ***//16X,'*** INTENSITY DATA ARE TO BE RE-', & 'NORMALIZED ***'//2X,'*** IT IS NOT RECOMMENDED TO SOLVE THE', & ' STRUCTURE BY PATTERSON ANALYSIS ***'/) RETURN END C ----------------------- SUBROUTINE COMMUL(I1,I2,I3,ICOM) C FIND A LEAST COMMON MULTIPLE FOR THE INCREMENTS C OF DIFFERENT INDICES DIMENSION IN(3) IN(1)=I1 IN(2)=I2 IN(3)=I3 ICOM=1 50 DO 300 J=2,256 DO 100 I=1,3 IF (MOD(IN(I),J).EQ.0) GO TO 150 100 CONTINUE 300 CONTINUE GO TO 400 150 DO 250 I=1,3 IF (MOD(IN(I),J).EQ.0) IN(I)=IN(I)/J 250 CONTINUE ICOM=ICOM*J DO 350 I=1,3 IF (IN(I).NE.1) GO TO 50 350 CONTINUE 400 RETURN END C ------------------------------------ SUBROUTINE GCD(EDX,ICX,INDEX,SCL,IP,IP6,NRLN) C FIND THE GREATEST COMMON DIVISORS FOR THE INCREMENTS C OF H K AND L OF THE STRONG REFLECTIONS DIMENSION EDX(3,256),ICX(3),SED(256),IP(8,3,5),IP6(6,5) DO 600 I=1,3 SIGE=0.0 DO 400 J=1,256 SED(J)=0.0 SIGE=SIGE+EDX(I,J) 400 CONTINUE DO 500 J=2,256 IF (EDX(I,J).LT.0.01) GO TO 500 DO 450 K=1,256 IF (EDX(I,K).LT.0.01) GO TO 450 IF (MOD(K,J).EQ.0) GO TO 450 SED(J)=SED(J)+EDX(I,K) 450 CONTINUE 500 CONTINUE SED(1) = 10000000.0 DO 550 J=2,256 IF (NRLN.EQ.0) GO TO 512 DO 510 N=1,NRLN IF (IP6(N,I).EQ.1.AND.IP6(N,4).EQ.J) EDX(I,J)=0.0 510 CONTINUE 512 IF (EDX(I,J).LT.0.01) GO TO 550 IF (SED(J)/SIGE.LT.SCL*0.6) GO TO 520 EDX(I,J)=0.0 GO TO 550 520 IF (SED(J).GT.SED(1)) GO TO 550 SED(1) = SED(J) ICX(I) = J 550 CONTINUE DO 570 J=2,256 IF (EDX(I,J).LT.0.01) GO TO 570 IF (MOD(J,ICX(I)).EQ.0) ICX(I)=J 570 CONTINUE 600 CONTINUE INDEX=1 DO 700 I=1,3 IF (ICX(I).GT.1) INDEX=INDEX+1 700 CONTINUE RETURN END C ------------------ SUBROUTINE GCD3(IN,IC,JUMP) C FIND THE GREATEST COMMON DIVISOR FOR THREE INTEGERS DIMENSION IN(3) JUMP=0 DO 200 I=1,9 J=10-I DO 100 K=1,3 IF (MOD(IN(K),J).NE.0) GO TO 200 100 CONTINUE GO TO 300 200 CONTINUE 300 IF (J.EQ.1.OR.MOD(J,IC).NE.0) JUMP=1 RETURN END C ---------------------------------- SUBROUTINE RCHECK(IN,ICOM,NRLN,JUMP,NB,IP6) C CHECK INDEX RELATIONSHIPS COMMON/REFLXOUT/INX(3,30000),EX(30000) COMMON/RESCALING/ISC,IBGR,MG,MIG(8),SCL,BT(9),SC(9),IP(8,3,5), 1 SCAL(8) DIMENSION IN(3),IP6(6,5) EBOT=0.0 ETOP=0.0 JUMP=0 DO 100 I=1,NB EBOT=EBOT+EX(I) INSM=IN(1)*INX(1,I)+IN(2)*INX(2,I)+IN(3)*INX(3,I) IF (MOD(INSM,ICOM).NE.0) ETOP=ETOP+EX(I) 100 CONTINUE IF (ETOP/EBOT.GT.SCL*0.6) GO TO 400 NRLN=NRLN+1 DO 300 I=1,3 IP6(NRLN,I)=IN(I) 300 CONTINUE IP6(NRLN,4)=ICOM JUMP=1 IP6(NRLN,5)=INT(1000000.0*ETOP/EBOT) 400 RETURN END C -------- SUBROUTINE SCRWRT C REWRITE THE SCRATCH FILE COMMON/REFLXIN/LH(600),LK(600),LL(600),FO(600),ID(600),RHO(600), 1 SIG(600),XKP(600) COMMON/WILSON/FLGW(30,9),FLGD(30,9),AVR(30,9),DCV(50,9),SLOPE(9), 1 FLGK(9),DEL(9),KS(9),EW(600),ED(600),EDP(600) COMMON/XXX/P(6),CX(9),NREF CALL CCPDPN(7,'SCRA7.TM','SCRATCH','U',80,0) c OPEN(7,FILE='SCRA7.TM',FORM='UNFORMATTED',STATUS='UNKNOWN') NREF=0 REWIND 8 100 READ (8) LH,LK,LL,FO,ID,EW,ED,RHO,SIG WRITE (7) LH,LK,LL,FO,ID,EW,ED,RHO,SIG DO 200 I=1,600 IF (FO(I).LT.0.0) GO TO 300 200 CONTINUE GO TO 100 300 REWIND 8 REWIND 7 400 READ (7) LH,LK,LL,FO,SIG,SIG,SIG,SIG,SIG CALL FCAL DO 500 I=1,600 IF (FO(I).LT.0.0) GO TO 600 500 CONTINUE GO TO 400 600 close (7) c 600 CLOSE (7,STATUS='DELETE') RETURN END C----------------------------------------------------------------------- SUBROUTINE CURVK(ESQ,RHO,ED,IG) C K-CURVE INTERPOLATION COMMON/WILSON/FLGW(30,9),FLGD(30,9),AVR(30,9),DCV(50,9),SLOPE(9), 1 FLGK(9),DEL(9),KS(9) COMMON/XXX/P(6),CX(9),NREF,NB,RHOMAX,MM,EMAX,EN,MZ,ER,ISIM,RHOCUT 1 ,RHOMIN,RHOLOW EI=RHO*DEL(IG) I=EI C FOR SMALL RHO THE CURVE IS EXTRAPOLATED FROM THE FIRST POINT C USING THE LEAST SQUARES SLOPE IF (I.LT.KS(IG)) SK = DCV(KS(IG),IG)*EXP(SLOPE(IG)*(FLOAT(I)/ 1 DEL(IG)-RHO)) C INTERPOLATION IF (I.GE.KS(IG).AND.I.LT.50) SK=DCV(I,IG)+(EI-FLOAT(I))*(DCV(I+1, 1 IG)-DCV(I,IG)) C FOR LARGE RHO THE CURVE IS EXTRAPOLATED FROM THE LAST POINT IF (I.GE.50) SK=DCV(50,IG)*EXP(SLOPE(IG)*(AVR(NB,IG)-RHO)) ESQ=ED/SK RETURN END C----------------------------------------------------------------------- SUBROUTINE DATAIN(MH,MK,ML) C READ REFLEXIONS FROM DATA FILE COMMON/REFLXIN/LH(600),LK(600),LL(600),FO(600),ID(600),RHO(600), 1 SIG(600),XKP(600) COMMON/REFLXOUT/IX(30000),EX(30000),FX(30000),SCMK(30000) COMMON/RESCALING/ISC,IBGR,MG,MIG(8),SCL,BT(9),SC(9),IP(8,3,5), 1 SCAL(8) COMMON/SYMMETRY/IS(3,3,24),TS(3,24),NSYM,PTS,KSYS,ICENT,LATT, 1 IAVE,EXTI COMMON/XXX/ P(6),CX(9),NREF,NB,RHOMAX,MM,EMAX,EN,MZ,ER,ISIM,RHOCUT REAL RSYM(4,4,64),CELLP(6) DIMENSION ADATA1(5),I1(3),I2(3),KL3(24) integer IFLAG(24),JFLAG(24) CHARACTER ITLE INTEGER IPOINT1(5) CHARACTER*30 LAB1(5),LSUSRJ(5) CHARACTER*1 CTYP1(5) LOGICAL EOF,LOGMSS(5) CHARACTER*4 KEY, CCVAL(20) CHARACTER*400 LINE CHARACTER NCH(7), LTYPEX*1,PGNAMX*10,SPGRNX*10 INTEGER NTOK, IBEG(20), IEND(20), ITYP(20), IDEC(20) REAL FVAL(20) DATA ADATA1/5*0.0/ DATA LAB1/'H','K','L','DF','SIGDF'/ DATA CTYP1/'H','H','H','D','Q'/ DATA IPOINT1/-1,-1,-1,2*0/ DATA NCH/'A','B','C','I','R','F','P'/ rewind (10) c CALL CCPDPN(10,'KEYWORD.TM','UNKNOWN','F',80,0) c OPEN(10,FILE='KEYWORD.TM',FORM='FORMATTED',STATUS='UNKNOWN') 10 READ(10,'(A132)',END=15) LINE NTOK = 20 CALL PARSER (KEY,LINE,IBEG,IEND,ITYP,FVAL,CCVAL, + IDEC,NTOK,EOF,.FALSE.) CALL CCPUPC(KEY) IF (KEY.NE.'LABI') GO TO 10 GOTO 20 15 CALL CCPERR(1,' ---- LABIN card missing ---- ') 20 NL = 5 CALL MTZINI ITOK = 2 CALL LKYSET(LAB1,NL,LSUSRJ,IPOINT1,ITOK,NTOK,LINE, + IBEG,IEND) CALL LKYIN(1,LAB1,NL,NTOK,LINE,IBEG,IEND) CALL LROPEN (1,'HKLIN',1,IFAIL) IF (IFAIL.EQ.1) CALL CCPERR(1,'Error opening HKLIN') if (latt.eq.0) then CALL LRSYMM(1,NSYM,RSYM) call lrsymi(1,nsympx,ltypex,nspgrx,spgrnx,pgnamx) DO 22 I=1,7 IF (LTYPEX.EQ.NCH(I)) IN1=I 22 CONTINUE LATT=MOD(IN1,7)+1 IF (LATT.LE.5) PTS=MIN0(2,LATT) IF (LATT.GE.6) PTS=LATT-3 IF (IN1.EQ.6) LATT=6 IF (IN1.EQ.5) LATT=7 nsym = nsym/pts do n = 1, nsym do i = 1, 3 ts(i,n) = rsym(i,4,n) do j = 1, 3 is(i,j,n) = nint(rsym(j,i,n)) end do end do end do ksys = 1 if (nspgrx.ge.3.and.nspgrx.le.15) ksys = 2 if (nspgrx.ge.16.and.nspgrx.le.74) ksys = 3 if (nspgrx.ge.75.and.nspgrx.le.142) ksys = 4 if (nspgrx.ge.143.and.nspgrx.le.167) ksys = 5 if (nspgrx.ge.168.and.nspgrx.le.194) ksys = 6 if (nspgrx.ge.195) ksys = 8 end if if (cx(1).eq.0.0) then CALL LRCELL(1,CELLP) do i = 1, 6 cx(i) = cellp(i) end do call incell end if CALL LRASSN (1,LAB1,NL,IPOINT1,CTYP1) K=0 100 CALL LRREFF(1,RESOL,ADATA1,EOF) IF (EOF) GOTO 450 CALL LRREFM(1,LOGMSS) IF (LOGMSS(4)) GOTO 100 IF (adata1(4).EQ.0.) GOTO 100 c IF (abs(adata1(4)).lt.(2.*adata1(5))) GOTO 100 KH=ADATA1(1) KK=ADATA1(2) KL=ADATA1(3) F=ABS(ADATA1(4)) SD=ADATA1(5) C GENERATE SYMMETRY RELATED REFLEXIONS AND FIND STANDARD ONE. MAXI=1 I1(1)=KH I1(2)=KK I1(3)=KL KCEN = 0 DO 250 J=1,NSYM IFLAG(J) = 1 JFLAG(J) = 1 DO 220 I0=1,3 I2(I0)=IS(I0,1,J)*I1(1) + IS(I0,2,J)*I1(2) + IS(I0,3,J)*I1(3) 220 CONTINUE IF((I2(1)+I1(1)).EQ.0.AND.(I2(2)+I1(2)).EQ.0.AND. + (I2(3)+I1(3)).EQ.0) KCEN=1 IND=65536*I2(1)+256*I2(2)+I2(3) IF (IND.LT.0) IFLAG(J) = -1 c IF (IND.LT.0.AND.IDTYPE.EQ.-1) JFLAG(J) = -1 KL3(J)=32896+IABS(IND) IF(J.EQ.1) GOTO 250 JM1=J-1 DO 240 I0=1,JM1 IF(KL3(I0).EQ.KL3(J)) KL3(J)=0 240 CONTINUE IF(KL3(J).GT.KL3(MAXI)) MAXI=J 250 CONTINUE C UNPACKING STANDARD REFLEXIONS IND=KL3(MAXI) KH=IND/65536 IF(IND.LT.0) KH=KH-1 IND=IND-65536*KH KK=IND/256 KL=IND-256*KK-128 KK=KK-128 K=K+1 LH(K)=KH IF (MH.LT.KH) MH=KH LK(K)=KK IF (MK.LT.KK) MK=KK LL(K)=KL IF (ML.LT.KL) ML=KL FO(K)=F SIG(K)=1. 300 CONTINUE IF(K.NE.600) GOTO 100 C CALCULATE RHO, EPSILON, MULTIPLICITY, SCATTERING FACTOR AND C CREATE SCRATCH FILE CALL FCAL K=0 GOTO 100 450 FO(K+1) = -1.0 CALL FCAL CALL LRCLOS (1) CLOSE (10) RETURN END C--------------------------------------------------------------- SUBROUTINE AVERAGE(FM,NPC) CHARACTER FM*80 COMMON/XXX/P(6) DIMENSION LHKL(4,30000),FODS(2,30000) DIMENSION KH(10),KK(10),KL(10),F(10),KD(10) C INPUT REFLECTION DATA K=0 80 READ (3,FM) (KH(I),KK(I),KL(I),F(I),KD(I),I=1,NPC) DO 100 I=1,NPC IF (F(I).LT.0.0) GO TO 120 K=K+1 LHKL(1,K)=KH(I) LHKL(2,K)=KK(I) LHKL(3,K)=KL(I) LHKL(4,K)=KD(I) FODS(1,K)=F(I) FODS(2,K)= * 0.5/SQRT(P(1)*FLOAT(KH(I)*KH(I))+P(2)*FLOAT(KK(I)*KK(I)) 1 +P(3)*FLOAT(KL(I)*KL(I))+P(4)*FLOAT(KH(I)*KK(I)) 2 +P(5)*FLOAT(KH(I)*KL(I))+P(6)*FLOAT(KK(I)*KL(I))) IF(K.EQ.30000) GO TO 120 100 CONTINUE GO TO 80 120 CONTINUE C DATA DEALING PROCESS CALL SORTD(6,LHKL,FODS,K) CALL DATAV(LHKL,FODS,K) NUMB=0 DO 180 LN=1,K IF (LHKL(4,LN).EQ.2) GO TO 180 NUMB=NUMB+1 LHKL(1,NUMB)=LHKL(1,LN) LHKL(2,NUMB)=LHKL(2,LN) LHKL(3,NUMB)=LHKL(3,LN) LHKL(4,NUMB)=LHKL(4,LN) FODS(1,NUMB)=FODS(1,LN) 180 CONTINUE CALL SORTD(1,LHKL,FODS,NUMB) NUMB=NUMB+1 LHKL(1,NUMB)=10 LHKL(2,NUMB)=10 LHKL(3,NUMB)=10 LHKL(4,NUMB)=1 FODS(1,NUMB)=-100.0 REWIND 3 WRITE(3,'(A80)') FM WRITE(3,FM) (LHKL(1,LN),LHKL(2,LN),LHKL(3,LN),FODS(1,LN), * LHKL(4,LN),LN=1,NUMB) REWIND 3 READ (3,'(A80)') FM RETURN END C ----------- SUBROUTINE DATAV(LHKL,FODS,NUMB) C SELECTING STANDARD REFLECTION DATA COMMON/SYMMETRY/IS(3,3,24),TS(3,24),NSYM,PTS,KSYS,ICENT,LATT, 1 IAVE,EXTI DIMENSION LHKL(4,600000),FODS(2,600000),MKS(3,3,48), * LLH(48),LLK(48),LLL(48) MP=NSYM*2 DO 40 I=1,NSYM DO 35 J=1,3 DO 30 M=1,3 MKS(M,J,I)=IS(M,J,I) MKS(M,J,I+NSYM)=-IS(M,J,I) 30 CONTINUE 35 CONTINUE 40 CONTINUE LN=1 100 KC=1 IF (LHKL(4,LN).EQ.2) GO TO 160 DO 110 I=1,MP LLH(I)=MKS(1,1,I)*LHKL(1,LN)+MKS(1,2,I)*LHKL(2,LN)+ * MKS(1,3,I)*LHKL(3,LN) LLK(I)=MKS(2,1,I)*LHKL(1,LN)+MKS(2,2,I)*LHKL(2,LN)+ * MKS(2,3,I)*LHKL(3,LN) LLL(I)=MKS(3,1,I)*LHKL(1,LN)+MKS(3,2,I)*LHKL(2,LN)+ * MKS(3,3,I)*LHKL(3,LN) 110 CONTINUE LH0=LLH(1) LK0=LLK(1) LL0=LLL(1) DO 120 I=1,MP-1 IF (LLH(I+1).LT.LH0) GO TO 120 IF (LLH(I+1).EQ.LH0.AND.LLK(I+1).LT.LK0) GO TO 120 IF (LLH(I+1).EQ.LH0.AND.LLK(I+1).EQ.LK0.AND.LLL(I+1).LT.LL0) * GO TO 120 LH0=LLH(I+1) LK0=LLK(I+1) LL0=LLL(I+1) 120 CONTINUE ICH1=0 ICH2=0 DO 140 J=LN+1,NUMB DJN=ABS(FODS(2,J)-FODS(2,LN)) IF (DJN.GT.1E-6) GO TO 150 ICH1=ICH1+1 DO 130 I=1,MP IF(LHKL(1,J).NE.LLH(I).OR.LHKL(2,J).NE.LLK(I).OR. * LHKL(3,J).NE.LLL(I)) GO TO 130 ICH2=ICH2+1 FODS(1,LN)=FODS(1,LN)+FODS(1,J) KC=KC+1 LHKL(4,J)=2 GO TO 140 130 CONTINUE 140 CONTINUE 150 FODS(1,LN)=FODS(1,LN)/KC LHKL(1,LN)=LH0 LHKL(2,LN)=LK0 LHKL(3,LN)=LL0 160 LN=LN+1 IF(LN.LT.NUMB) GOTO 100 RETURN END C------------------------------------------------------- SUBROUTINE SORTD(IFLAG,LHKL,FODS,NUMB) DIMENSION LHKL(4,600000),FODS(2,600000),LH(4),FO(2) INT=NUMB 100 INT=INT/2 IF(2*(INT/2).EQ.INT) INT=INT-1 IFIN=NUMB-INT DO 1000 II=1,IFIN I=II J=I+INT IF (IFLAG.LE.4) GO TO 150 IF (FODS(2,I).LE.FODS(2,J)) GO TO 1000 A1= FODS(2,J) GO TO 160 150 IF (LHKL(1,I).LT.LHKL(1,J)) GO TO 1000 IA1= LHKL(1,J) IF (LHKL(1,I).EQ.LHKL(1,J).AND.LHKL(2,I).LT.LHKL(2,J)) GO TO 1000 IA2= LHKL(2,J) IF (LHKL(1,I).EQ.LHKL(1,J).AND.LHKL(2,I).EQ.LHKL(2,J).AND. * LHKL(3,I).LT.LHKL(3,J)) GO TO 1000 IA3= LHKL(3,J) 160 DO 200 K=1,4 LH(K)=LHKL(K,J) 200 CONTINUE DO 300 K=1,2 FO(K)=FODS(K,J) 300 CONTINUE 400 DO 500 K=1,4 LHKL(K,J)=LHKL(K,I) 500 CONTINUE DO 600 K=1,2 FODS(K,J)=FODS(K,I) 600 CONTINUE J=I I=I-INT IF(I.LE.0) GOTO 700 IF (IFLAG.LE.4) GO TO 650 IF(FODS(2,I).GT.A1) GOTO 400 GO TO 700 650 IF(LHKL(1,I).GT.IA1) GOTO 400 IF(LHKL(1,I).EQ.IA1.AND.LHKL(2,I).GT.IA2) GO TO 400 IF(LHKL(1,I).EQ.IA1.AND.LHKL(2,I).EQ.IA2.AND. * LHKL(3,I).GT.IA3) GO TO 400 700 DO 800 K=1,4 LHKL(K,J)=LH(K) 800 CONTINUE DO 900 K=1,2 FODS(K,J)=FO(K) 900 CONTINUE 1000 CONTINUE IF(INT.GT.1) GOTO 100 RETURN END C ----------------------------------------------------------------- C THIS PROGRAM IS USED FOR CHECKING AND REJECTING THE REFLEXIONS C WITH SYSTEMATIC EXTINTION SUBROUTINE REJECT(IEXT,I1) COMMON/SYMMETRY/IS(3,3,24),TS(3,24),NSYM,PTS,KSYS,ICENT,LATT, 1 IAVE,EXTI DIMENSION I1(3),T(3) IHKL=I1(1)+I1(2) C LATT. CENTERED (P:500 A:400 B:200 C:450 I:300 F:350 R:100) GO TO (500,400,200,450,300,350,100) LATT 100 IF (MOD((IHKL+I1(1)+I1(3)),3).EQ.0) GO TO 500 GO TO 1000 200 IHKL=I1(1)+I1(3) GO TO 450 300 IHKL=IHKL+I1(3) GO TO 450 350 IF (MOD(IHKL,2).EQ.0) GO TO 400 GO TO 1000 400 IHKL=I1(2)+I1(3) 450 IF (MOD(IHKL,2).EQ.0) GO TO 500 GO TO 1000 500 IF (EXTI.LT.0.01) GO TO 1100 TEMP1=-999999.0 RECY=0.0 DO 980 J=1,NSYM DO 900 L=1,3 T(L)=IS(L,1,J)*I1(1) + IS(L,2,J)*I1(2) + IS(L,3,J)*I1(3) 900 CONTINUE IF (49.5-AMAX1(ABS(T(1)),ABS(T(2)))) 1100,1100,950 950 TEMP2=T(1)+100.00*(T(2)+100.0*T(3)) IF (ICENT.EQ.1) TEMP2=ABS(TEMP2) IF (TEMP2-TEMP1+0.01) 980,960,960 960 ANG=6.283185*(I1(1)*TS(1,J)+I1(2)*TS(2,J)+I1(3)*TS(3,J)) SUM=COS(ANG) C FOR THE ACENTRIC CASE IF (ICENT.EQ.0) SUM=SUM+100.0*SIN(ANG) ANG=SUM RECY=RECY+ANG IF (TEMP2-TEMP1-0.05) 980,970,970 970 RECY=ANG TEMP1=TEMP2 980 CONTINUE IF (ABS(RECY)-0.01) 1000,1100,1100 C FOR THE CASE OF SYSTEMATIC EXTINCTION 1000 IEXT=1 RETURN C FOR THE CASE OF NON-SYSTEMATIC EXTINCTION 1100 IEXT=0 RETURN END C----------------------------------------------------------------------- SUBROUTINE DEBYE(NMEM,NTOT,CR,POP) C CALCULATE SCATTERING FACTORS FOR RANDOM GROUPS (NINF = 2) COMMON/ATMFACTOR/AL(8),AS(8),BL(8),BS(8),CL(8),CS(8),DL(8),DS(8), 1 EL(8),NW(8),NO(8),NK,NAT,F(9) COMMON/ATMGROUP/NINF(10),NAG(10),X(5000),Y(5000),Z(5000),NZ(5000), 1 Q(5000),U11(5000),U22(5000),U33(5000),U23(5000), 2 U13(5000),U12(5000) COMMON/SCATFACTOR/GIS(142),GIW(142),NGP,NDIFF,ISOL,X0(3) COMMON/UNIT1/ITLE(80),LIST,PI CHARACTER ITLE DIMENSION G(142),CR(9) FP=4.0*PI IF (CR(1).LT.0.01) GO TO 15 WRITE (6,5) (CR(I),I=1,6) 5 FORMAT(36H GIVEN COORDINATES REFER TO THE CELL/2X,4H A =,F6.2, 1 2X,3HB =,F6.2,2X,3HC =,F6.2,2X,7HALPHA =,F6.1,2X,6HBETA =, 2 F6.1,2X,7HGAMMA =,F6.1) C CONVERT TO ORTHOGONAL COORDINATES CALL VOL(CR,V) R1=CR(3)*(CR(4)-CR(5)*CR(6))/CR(9) R2=CR(3)*V/CR(9) DO 10 I=NMEM,NTOT X(I)=CR(1)*X(I)+CR(2)*CR(6)*Y(I)+CR(3)*CR(5)*Z(I) Y(I)=CR(2)*CR(9)*Y(I)+R1*Z(I) Z(I)=R2*Z(I) 10 CONTINUE 15 D=1.0 DO 100 I=1,142 G(I)=0.0 T=0.01*FLOAT(I-1) TT=T*T DO 20 J=1,NK F(J)=AL(J)*EXP(-AS(J)*TT)+BL(J)*EXP(-BS(J)*TT)+CL(J)*EXP(-CS(J) 1 *TT)+DL(J)*EXP(-DS(J)*TT)+EL(J) 20 CONTINUE C SUM OVER INTERATOMIC VECTORS DO 50 K=NMEM,NTOT KT=NZ(K) NT=K+1 G(I)=G(I)+F(KT)*F(KT) IF (K.EQ.NTOT) GO TO 50 DO 40 L=NT,NTOT LT=NZ(L) IF (I.EQ.1) GO TO 30 D=FP*T*SQRT((X(K)-X(L))**2+(Y(K)-Y(L))**2+(Z(K)-Z(L))**2) D=SIN(D)/D 30 G(I)=G(I)+2.0*F(KT)*F(LT)*D 40 CONTINUE 50 CONTINUE G(I)=SQRT(G(I)) 100 CONTINUE WRITE (6,120) (G(I),I=1,101,2) 120 FORMAT(/1X,26X,'GROUP SCATTERING FACTORS AT INTERVALS OF 0.02 IN', 1 'SIN(THETA)/LAMBDA'/(1H ,17F7.2)) DO 140 I=1,142 GIS(I)=GIS(I)+G(I)*G(I)*POP 140 CONTINUE RETURN END ************************************************************************ * * * EEEEEEE CCCCC A L * * E C C A A L * * E C A A L * * EEEEE C A A L * * E C AAAAAAA L * * E C C A A L * * EEEEEEE CCCCC A A LLLLLLL * * * * CALCULATE FINAL E-VALUES AND RESCALED F'S * * OUTPUT REFLEXIONS FOR PHASE, PREPARE TABLES OF STATISTICS * * VERSION 1998 * ************************************************************************ SUBROUTINE ECAL COMMON/ATMGROUP/NINF(10),NAG(10),X(5000),Y(5000),Z(5000),NZ(5000), 1 Q(5000),U11(5000),U22(5000),U33(5000),U23(5000), 2 U13(5000),U12(5000) COMMON/REFLXIN/LH(600),LK(600),LL(600),FO(600),ID(600),RHO(600), 1 SIG(600),XKP(600) COMMON/REFLXOUT/IX(30000),EX(30000),FX(30000),SCMK(30000) COMMON/RESCALING/ISC,IBGR,MG,MIG(8),SCL,BT(9),SC(9),IP(8,3,5), 1 SCAL(8),RSTT COMMON/SCATFACTOR/GIS(142),GIW(142),NGP,NDIFF,ISOL,X0(3) COMMON/STATISTICS/VST(10,5),NST(5),ZT(25,5),EE(10),MULT,IND,NZR, 1 TMUL COMMON/SYMMETRY/IS(3,3,24),TS(3,24),NSYM,PTS,KSYS COMMON/UNIT1/ITLE(80),LIST,PI,KCURV COMMON/WILSON/FLGW(30,9),FLGD(30,9),AVR(30,9),DCV(50,9),SLOPE(9), 1 FLGK(9),DEL(9),KS(9),EW(600),ED(600),EDP(600) COMMON/XXX/P(6),CX(9),NREF,NB,RHOMAX,MM,EMAX,EN,MZ,ER,ISIM,RHOCUT, 1 RHOMIN,RHOLOW,ENG(8),ERG(8),KNOWN,IPATH,VOLUME DIMENSION RHR(10),NHR(10),NU(25),STL(10),INP(600),NOB(600), 1 INPA(600),IXW(30000),EXW(30000),FXW(30000),EH(600) C TABLES OF THEORETICAL DISTRIBUTIONS CHARACTER ITLE,ICHR(80),IOBS,IUNOBS,NOB DIMENSION AVA(10),AVC(10),AVH(10) C DIMENSION CPH(25) DATA AVA/0.886, 1.0,1.329, 2.0,3.323, 6.0,0.736,1.0,2.0,2.415/ DATA AVC/0.798, 1.0,1.596, 3.0,6.383,15.0,0.968,2.0,8.0,8.691/ DATA AVH/0.718, 1.0,1.916,4.5,12.26,37.5,1.145,3.5,26.0,26.903/ C DATA CPH/0.368,0.463,0.526,0.574,0.612,0.643,0.670,0.694,0.715, C 1 0.733,0.765,0.791,0.813,0.832,0.848,0.863,0.875,0.886,0.896, C 2 0.905,0.913,0.920,0.926,0.932,0.938/ DATA IOBS,IUNOBS/1H ,1HU/ DATA ICHR/19*1H ,1HA,1HL,1HL,1H ,1HD,1HA,1HT,1HA,4*1H ,1HH,1HK, 1 1HL,4*1H ,1H0,1HK,1HL,17*1H ,1HA,1HC,1HE,1HN,1HT,1H ,1H ,1HC, 2 1HE,1HN,1HT,1H ,1H ,1HH,1H-,1HC,1HE,1HN,1HT,3*1H / C SET INITIAL VALUES IF (ISOL.EQ.1) MMS=MM IF (ISOL.EQ.1) MM=15000 DO 20 I=1,10 RHR(I)=0.0 NHR(I)=0 DO 10 J=1,5 VST(I,J)=0.0 10 CONTINUE 20 CONTINUE DO 40 I=1,5 NST(I)=0 DO 30 J=1,25 ZT(J,I)=0.0 30 CONTINUE 40 CONTINUE DO 50 I=1,25 NU(I)=0 50 CONTINUE NRW=0 LIM=30000 NC=0 NS=0 NL=0 IF (RSTT.GT.1.0) GO TO 54 IF (RSTT.GE.0.0) GO TO 53 MGW=0 DO 52 I=1,MG IF (MIG(I).EQ.-1) MGW=MGW+1 52 CONTINUE IF (MGW.EQ.0) GO TO 54 RSTT=FLOAT(MG-MGW)/FLOAT(MG) 53 LIMW=30000 NCW=0 NSW=0 NLW=0 ENW=EN ERW=ER MMW=INT(FLOAT(MM)*(1.0-RSTT)) MM=MM-MMW MZW=INT(FLOAT(MZ)*(1.0-RSTT)) MZ=MZ-MZW 54 KG=1 SCF=SQRT(SC(1)) RR=10.0/SQRT(RHOMAX) IF (LIST.EQ.1) WRITE (6,56) NREF 56 FORMAT(//1X,7HLIST OF,I5,1X,'REFLEXIONS.'/' F IS ON AN ABSOLUTE ', 1 'SCALE AND THE E-VALUES ARE SUITABLE FOR INPUT TO PHASE'/ 2 1X,1X,3(9H H K L,5X,1HF,4X,1HE,4X)) REWIND 8 60 READ (8) LH,LK,LL,FO,ID,EW,ED,RHO,SIG,EDP DO 250 I=1,600 IG=MOD(ID(I),100) IF (IBGR.EQ.1) KG=IG INPA(I)=32896 MF=I C TEST FOR END OF DATA IF (FO(I).LT.0.0) GO TO 300 D=EXP(BT(IG)*RHO(I)) D=D*SCAL(IG) MULT=ID(I)/10000 IE=(ID(I)-10000*MULT)/100 IF (KCURV.NE.1) ESQ=D*ED(I) IF (KCURV.EQ.1) CALL CURVK(ESQ,RHO(I),ED(I),KG) C CALCULATE E AND RESCALED F ED(I)=SQRT(ESQ) FO(I)=FO(I)*SCF C PACKING INPA(I)=65536*LH(I)+256*(LK(I)+128)+LL(I)+128 INP(I)=INPA(I)*32+ISIGN(1,INPA(I))*IG C SET OBS/UNOBS FLAG NOB(I)=IOBS IF (SIG(I).LT.0.0) NOB(I)=IUNOBS IF (SIG(I).GE.0.0 .AND. ED(I).LT.0.01) GO TO 110 IF (RHO(I).GT.RHOCUT) GO TO 110 IF (RHO(I).LT.RHOLOW) GO TO 110 C STORE REFLEXIONS FOR PHASE IF (RSTT.GT.1.0.OR.MIG(IG).GT.0) GO TO 85 IF (SIG(I).LT.0.0) GO TO 70 NLW=NLW+1 IF (ED(I).GT.ENW.AND.ED(I).LT.EMAX) GO TO 80 NLW=NLW-1 70 ED(I)=ED(I)+0.2*RHO(I) IF (ED(I).GT.ERW) GO TO 110 NSW=NSW+1 80 NCW=NCW+1 IXW(NCW)=INP(I) EXW(NCW)=ED(I) FXW(NCW)=FO(I) IF (NCW.NE.LIMW) GO TO 110 CALL SORT(EXW,FXW,IXW,30000) IF (MZW.GT.0) ERW=AMIN1(ERW,EXW(30001-MZW)) ENW=AMAX1(ENW,EXW(MMW)) NCW=MIN0(MMW,NLW) LIMW=30000-MZW GO TO 110 85 IF (SIG(I).LT.0.0) GO TO 90 NL=NL+1 IF (ED(I).GT.EN.AND.ED(I).LT.EMAX) GO TO 100 NL=NL-1 90 ED(I)=ED(I)+0.2*RHO(I) IF (ED(I).GT.ER) GO TO 110 NS=NS+1 100 NC=NC+1 IX(NC)=INP(I) EX(NC)=ED(I) FX(NC)=FO(I) IF (NC.NE.LIM) GO TO 110 CALL SORT(EX,FX,IX,30000) IF (MZ.GT.0) ER=AMIN1(ER,EX(30001-MZ)) EN=AMAX1(EN,EX(MM)) NC=MIN0(MM,NL) LIM=30000-MZ C WORK OUT FINAL STATISTICS 110 TMUL=FLOAT(MULT) C DISTRIBUTION OF E WITH SIN(THETA)/LAMBDA N=MIN0(10,INT(1.0+RR*SQRT(RHO(I)))) NHR(N)=NHR(N)+MULT RHR(N)=RHR(N)+ESQ*TMUL NZR=INT(10.0*ESQ)+1 IF (NZR.GT.10) NZR=10+(NZR-9)/2 EE(1)=ED(I) DO 120 J=2,6 EE(J)=EE(J-1)*ED(I) 120 CONTINUE EE(7)=ESQ-1.0 EE(8)=EE(7)*EE(7) EE(9)=EE(8)*EE(7) EE(10)=ABS(EE(9)) EE(7)=ABS(EE(7)) DO 130 J=1,10 EE(J)=TMUL*EE(J) 130 CONTINUE C ADD FUNCTIONS OF E TO APPROPRIATE ZONES IND=1 J=LH(I) K=LK(I) L=LL(I) CALL ADD(1) GO TO (140,140,140,150,160,160,170,180),KSYS C TRICLINIC, MONOCLINIC AND ORTHORHOMBIC 140 IF (J.EQ.0) CALL ADD(3) IF (K.EQ.0) CALL ADD(4) IF (L.EQ.0) CALL ADD(5) GO TO 200 C TETRAGONAL 150 IF (J.EQ.0.OR.K.EQ.0) CALL ADD(3) IF (IABS(J).EQ.IABS(K)) CALL ADD(4) IF (L.EQ.0) CALL ADD(5) GO TO 200 C TRIGONAL, HEXAGONAL AND RHOMBOHEDRAL INDEXED ON HEXAGONAL AXES 160 IF (J.EQ.0.OR.K.EQ.0.OR.J+K.EQ.0) CALL ADD(3) IF (J.EQ.K.OR.J+2*K.EQ.0.OR.2*J+K.EQ.0) CALL ADD(4) IF (L.EQ.0) CALL ADD(5) GO TO 200 C PRIMITIVE RHOMBOHEDRAL 170 IF (J.EQ.K.OR.J.EQ.L.OR.K.EQ.L) CALL ADD(3) IF (L.EQ.2*K-J.OR.K.EQ.2*J-L.OR.J.EQ.2*L-K) CALL ADD(4) IF (J+K+L.EQ.0) CALL ADD(5) GO TO 200 C CUBIC 180 IF (J.EQ.0.OR.K.EQ.0.OR.L.EQ.0) CALL ADD(3) IF (IABS(J).EQ.IABS(K).OR.IABS(J).EQ.IABS(L).OR. 1 IABS(K).EQ.IABS(L)) CALL ADD(4) C H,H,2H IS IN TWO PRINCIPAL ZONES BUT NOT ON A PRINCIPAL AXIS IF (IND.EQ.4) IND=0 IF (IABS(L).EQ.IABS(J+K).OR.IABS(K).EQ.IABS(J+L).OR. 1 IABS(J).EQ.IABS(K+L)) CALL ADD(5) C REFLEXIONS NOT BELONGING TO PRINCIPAL ZONES 200 IF (IND.EQ.1) CALL ADD(2) C DISTRIBUTION OF E FOR COMPLETE DATA NET=MIN0(25,INT(10.0*ED(I))) IF (NET.EQ.0) GO TO 220 NU(NET)=NU(NET)+1 220 NRW=NRW+MULT 250 CONTINUE 300 IF (FO(MF).LT.0.0) MF=MF-1 C LIST REFLEXIONS IF REQUIRED IF (LIST.EQ.1.AND.MF.NE.0) WRITE (6,306) (LH(K),LK(K),LL(K),FO(K), 1 ED(K),NOB(K),K=1,MF) C FORMAT FOR REFLEXION LIST - A WIDER LINEPRINTER MAY ALLOW MORE C REFLEXIONS TO BE OUTPUT ON ONE LINE - THUS SAVING PAPER 306 FORMAT(3(2H ,3I3,F7.1,F5.2,1X,A1)) C TEST FOR END OF DATA IF (MF.EQ.600) GO TO 60 C OUTPUT STATISTICS RR=1.0/RR DO 320 I=1,10 STL(I)=RR*FLOAT(I) IF (NHR(I).GT.0) RHR(I)=RHR(I)/FLOAT(NHR(I)) DO 310 J=1,5 IF (NST(J).GT.0) VST(I,J)=VST(I,J)/FLOAT(NST(J)) 310 CONTINUE 320 CONTINUE WRITE (6,330) STL,RHR,NHR 330 FORMAT(///1X,32X,16HFINAL STATISTICS//1H ,17X, 1 43HDISTRIBUTION OF E**2 WITH SIN(THETA)/LAMBDA/8H SIN/LAM, 2 F6.4,9F7.4/1H ,1X,4HE**2,2X,F6.4,9F7.4/1H ,6HNUMBER,10I7) WRITE (6,340) 340 FORMAT(//1X,32X,14HAVERAGE VALUES, 1 //1X,4X,7HAVERAGE,19X,'EXPERIMENTAL',19X,'THEORETICAL') IF (KSYS.LE.6) THEN ICHR(53)='H' ICHR(54)='K' ICHR(55)='0' IF (KSYS.LE.3) THEN ICHR(46)='H' ICHR(47)='0' ICHR(48)='L' ELSE ICHR(46)='H' ICHR(47)='H' ICHR(48)='L' ENDIF ELSE ICHR(50)='H' ICHR(51)=',' ICHR(52)='K' ICHR(53)=',' ICHR(54)='-' ICHR(55)='H' ICHR(56)='-' ICHR(57)='K' IF (KSYS.EQ.7) THEN ICHR(35)='H' ICHR(36)=',' ICHR(37)='K' ICHR(38)='2' ICHR(39)='K' ICHR(40)='-' ICHR(41)='H' ENDIF ENDIF WRITE (6,420) ICHR 420 FORMAT(80A1/) WRITE (6,430) ((VST(I,J),J=1,5),AVA(I),AVC(I),AVH(I),I=1,10) 430 FORMAT(1H ,5X,6HMOD(E),8X,8F7.3/ 1 1H ,6X,4HE**2,9X,8F7.3/ 2 1H ,6X,4HE**3,9X,8F7.3/ 3 1H ,6X,4HE**4,9X,8F7.3/ 4 1H ,6X,4HE**5,9X,8F7.3/ 5 1H ,6X,4HE**6,9X,8F7.3/ 6 1H ,2X,11HMOD(E**2-1),6X,8F7.3/ 7 1H ,2X,11H(E**2-1)**2,6X,8F7.3/ 8 1H ,2X,11H(E**2-1)**3,6X,8F7.3/ 9 1H ,16H(MOD(E**2-1))**3,3X,8F7.3) WRITE (6,440) NST 440 FORMAT(1X,20HWEIGHTED SAMPLE SIZE,I6,4I7) WRITE (6,520) 520 FORMAT(//18X,'DISTRIBUTION OF E - NUMBER OF ',14HE'S .GT. LIMIT) DO 540 I=1,25 AVR(I,IG)=0.1*FLOAT(I) II=I+1 IF (II.GT.25) GO TO 540 DO 530 J=II,25 NU(I)=NU(I)+NU(J) 530 CONTINUE 540 CONTINUE WRITE (6,542) (AVR(I,IG),I=7,16),(NU(I),I=7,16), 1 (AVR(I,IG),I=17,25),(NU(I),I=17,25) 542 FORMAT(1H ,7X,3HE ,10F6.1/1H ,6X,4HNO. ,10I6// 1 1H ,7X,3HE ,9F6.1,6X/1H ,6X,4HNO. ,9I6) C OUTPUT REFLEXIONS FOR PHASE MR=30000 IF (LIM.EQ.30000) MR=NC CALL SORT(EX,FX,IX,MR) MM=MIN0(MM,NL) MZ=MIN0(MZ,NS) IF (RSTT.GT.1.0.OR.RSTT.LT.0.0) GO TO 562 MRW=30000 IF (LIMW.EQ.30000) MRW=NCW CALL SORT(EXW,FXW,IXW,MRW) MMW=MIN0(MMW,NLW) MZW=MIN0(MZW,NSW) I1=MR-MZ DO 544 I=1,MZ I2=MM+I I3=I1+I IX(I2)=IX(I3) EX(I2)=EX(I3) FX(I2)=FX(I3) 544 CONTINUE I1=MM+MZ DO 546 I=1,MMW I2=I1+I IX(I2)=IXW(I) EX(I2)=EXW(I) FX(I2)=FXW(I) 546 CONTINUE I1=I1+MMW I2=MRW-MZW DO 548 I=1,MZW I3=I1+I I4=I2+I IX(I3)=IXW(I4) EX(I3)=EXW(I4) FX(I3)=FXW(I4) 548 CONTINUE MM=MM+MMW MZ=MZ+MZW MR=MM+MZ CALL SORT(EX,FX,IX,MR) 562 IF (ISOL.NE.1) GO TO 570 MMS=MIN0(MMS,MM) NF=0 DO 565 IGP=1,NGP NNS=NF+1 NF=NF+NAG(IGP) IF (NINF(IGP).NE.8) GO TO 565 CALL SOLENA(NNS,NF,MMS) GO TO 570 565 CONTINUE 570 WRITE (6,572) MM 572 FORMAT(///1H ,I6,41H LARGEST E-VALUES WRITTEN TO OUTPUT FILE) IF (RSTT.GT.1.0.OR.RSTT.LT.0.0) GO TO 580 MM2=MM-MMW WRITE (6,576) MM2 576 FORMAT(8X,10H(INCLUDING,I4,15H SYSTEMATICALLY, 1 20H STRONG REFLECTIONS)) MM2=MM-MMW 580 DO 583 I=1,MG SCAL(I)=SCAL(I)/SC(1) 583 CONTINUE CALL OUTPUT(0,MM) MS=MR-MZ WRITE (6,600) MZ 600 FORMAT(/1X,I6,42H SMALLEST E-VALUES WRITTEN TO OUTPUT FILE) IF (RSTT.GT.1.0.OR.RSTT.LT.0.0) GO TO 620 MZ2=MZ-MZW c WRITE (6,574) MZ2 IF (RSTT.GT.1.0.OR.RSTT.LT.0.0) GO TO 620 MZ2=MZ-MZW 620 CALL OUTPUT(MS,MR) RETURN END C -------- SUBROUTINE ADD(N) C SUMS FOR REFLEXION IN ZONE N COMMON/STATISTICS/VST(10,5),NST(5),ZT(25,5),EE(10),MULT,IND,NZR, 1 TMUL IS=1 NT=N IF (N.LE.2.OR.IND.LE.1) GO TO 20 C REFLEXION IS ON PRINCIPAL AXIS - THEREFORE IGNORE IT NT=IND IS=-1 20 S=FLOAT(IS) DO 30 I=1,10 VST(I,NT)=VST(I,NT)+S*EE(I) 30 CONTINUE NST(NT)=NST(NT)+IS*MULT IF (NZR.LE.25) ZT(NZR,NT)=ZT(NZR,NT)+S*TMUL IND=N RETURN END C ---------------- SUBROUTINE SORT(A,B,IX,N) C SORT ON A DIMENSION A(N),B(N),IX(N) INT=2 10 INT=2*INT IF (INT.LT.N) GO TO 10 INT=MIN0(N,(3*INT)/4-1) 20 INT=INT/2 IFIN=N-INT DO 70 II=1,IFIN I=II J=I+INT IF (A(I).GE.A(J)) GO TO 70 T=A(J) X=B(J) L=IX(J) 40 A(J)=A(I) B(J)=B(I) IX(J)=IX(I) J=I I=I-INT IF (I.LE.0) GO TO 60 IF (A(I).LT.T) GO TO 40 60 A(J)=T B(J)=X IX(J)=L 70 CONTINUE IF (INT.GT.1) GO TO 20 RETURN END ************************************************************************ * * * FFFFFFF CCCCC A L * * F C C A A L * * F C A A L * * FFFFF C A A L * * F C AAAAAAA L * * F C C A A L * * F CCCCC A A LLLLLLL * * * * CALCULATE RHO, EPSILON, MULTIPLICITY AND SCATTERING FACTOR * * CREATE THE SCRATCH FILE (CH.=8), PREPARE FILE FOR WEIGHTED FOURIER * * VERSION 1998 * ************************************************************************ SUBROUTINE FCAL COMMON/ATMGROUP/NINF(10),NAG(10),X(5000),Y(5000),Z(5000),NZ(5000), 1 Q(5000),U11(5000),U22(5000),U33(5000),U23(5000), 2 U13(5000),U12(5000) COMMON/REFLXIN/LH(600),LK(600),LL(600),FO(600),ID(600),RHO(600), 1 FC(600),XKP(600) COMMON/REFLXOUT/IX(30000),EX(30000),FX(30000),SCMK(30000) COMMON/RESCALING/ISC,IBGR,MG,MIG(8),SCL,BT(9),SC(9),IP(8,3,5), 1 SCAL(8) COMMON/SCATFACTOR/GIS(142),GIW(142),NGP,NDIFF,ISOL,X0(3) COMMON/SYMMETRY/IS(3,3,24),TS(3,24),NSYM,PTS,KSYS,ICENT,LATT COMMON/UNIT1/ ITLE(80),LIST,PI,KCURV,IPAT,IAPA,MOV,MFP,NAU COMMON/WILSON/FLGW(30,9),FLGD(30,9),AVR(30,9),DCV(50,9),SLOPE(9), 1 FLGK(9),DEL(9),KS(9),EW(600),ED(600),EDP(600) COMMON/XXX/P(6),CX(9),NREF,NB,RHOMAX,MM,EMAX,EN,MZ,ER,ISIM,RHOCUT, 1 RHOMIN,RHOLOW,ZCG(8),ZOG(8),KNOWN,IPATH,VOLUME,ZC,ZO DIMENSION I1(3),I2(3),AC(600),BC(600) CHARACTER ITLE DO 300 I=1,600 IF (FO(I).LT.0.0) GO TO 310 NREF=NREF+1 RHO(I)=P(1)*FLOAT(LH(I)*LH(I))+P(2)*FLOAT(LK(I)*LK(I)) 1 +P(3)*FLOAT(LL(I)*LL(I))+P(4)*FLOAT(LH(I)*LK(I)) 2 +P(5)*FLOAT(LH(I)*LL(I))+P(6)*FLOAT(LK(I)*LL(I)) RHOMAX=AMAX1(RHOMAX,RHO(I)) RHOMIN=AMIN1(RHOMIN,RHO(I)) C COMPUTE EPSILON AND MULTIPLICITY BY GENERATING EQUIVALENT C REFLEXIONS C EPSILON = NUMBER OF TIMES SAME REFLEXION APPEARS IN LIST C MULTIPLICITY = NUMBER DIFFERENT REFLEXIONS IN LIST EPS=1.0 MULT=1 I1(1)=LH(I) I1(2)=LK(I) I1(3)=LL(I) C IN TRICLINIC SPACE GROUPS EPS = 1.0 AND MULT = 1 IF (NSYM.EQ.1) GO TO 60 K1=65536*I1(1)+256*(I1(2)+128)+I1(3)+128 IK1=65792-K1 DO 50 J=2,NSYM DO 40 L=1,3 I2(L)=IS(L,1,J)*I1(1) + IS(L,2,J)*I1(2) + IS(L,3,J)*I1(3) 40 CONTINUE K2=65536*I2(1)+256*(I2(2)+128)+I2(3)+128 IF (K2.EQ.K1) EPS=EPS+1.0 IF (ICENT.NE.0.AND.K2.EQ.IK1) EPS=EPS+1.0 IF (K2.EQ.K1.OR.K2.EQ.IK1) MULT=MULT+1 50 CONTINUE 60 FF=FO(I)*FO(I)/PTS ID(I)=NSYM/MULT IF (ISC.EQ.2) GO TO 92 C DETERMINE INDEX GROUP (FOR RESCALING) MG=8 IF (KSYS.GT.3) MG=6 IF (KSYS.GE.7) GO TO 90 LG=MOD(IABS(LL(I)),2) IF (KSYS.GE.5) GO TO 80 KG=MOD(IABS(LK(I)),2) JG=MOD(IABS(LH(I)),2) C TRICLINIC, MONCLINIC AND ORTHORHOMBIC IF (KSYS.LE.3) IG=JG+2*KG+4*LG C TETRAGONAL IF (KSYS.EQ.4) IG=JG+KG+3*LG GO TO 100 C TRIGONAL, HEXAGONAL AND RHOMBOHEDRAL INDEXED ON HEXAGONAL AXES 80 IG=3*LG IF (MOD(LH(I),3).EQ.0) IG=IG+1 IF (MOD(LK(I),3).EQ.0.OR.MOD(LH(I)+LK(I),3).EQ.0) IG=IG+1 GO TO 100 C CUBIC AND PRIMITIVE RHOMBOHEDRAL 90 IG=3*MOD(IABS(LH(I)+LK(I)+LL(I)),2) IF (MOD(LH(I)-LL(I),3).EQ.0) IG=IG+1 IF (MOD(LK(I)-LL(I),3).EQ.0.OR.MOD(LH(I)-LK(I),3).EQ.0) IG=IG+1 GO TO 100 C DETERMINE INDEX GROUP FOR SPECIAL RESCALING 92 DO 98 J=1,MG DO 96 L=1,3 K=IABS(IP(J,L,1)*LH(I)+IP(J,L,2)*LK(I)+IP(J,L,3)*LL(I) 1 -IP(J,L,5)) IF (IP(J,L,4).EQ.0) GO TO 94 IF (MOD(K,IP(J,L,4)).NE.0) GO TO 98 GO TO 96 94 IF (K.NE.0) GO TO 98 96 CONTINUE IG=J-1 GO TO 100 98 CONTINUE IG=MG C PACK SYMMETRY FUNCTIONS FOR LATER USE 100 IG=IG+1 ID(I)=10000*ID(I)+100*INT(EPS+0.5)+IG C LOOK UP SCATTERING FACTOR TABLES GENERATED BY ATMCOEF SINTH=100.0*SQRT(RHO(I)) IND=MAX0(2,INT(SINTH+1.5)) FRAC=SINTH-FLOAT(IND-1) BF=0.5*(GIW(IND+1)-GIW(IND-1)) AF=BF+GIW(IND-1)-GIW(IND) WFORM=AF*FRAC*FRAC+BF*FRAC+GIW(IND) C 'WILSON' STRUCTURE FACTOR EW(I)=FF/(WFORM*EPS) ED(I)=EW(I) IF (KCURV.EQ.-1) ED(I)=FF BF=0.5*(GIS(IND+1)-GIS(IND-1)) AF=BF+GIS(IND-1)-GIS(IND) FORM=AF*FRAC*FRAC+BF*FRAC+GIS(IND) IF (KNOWN.EQ.1.OR.NDIFF.EQ.1) GO TO 105 C 'DEBYE' STRUCTURE FACTOR ED(I)=FF/(FORM*EPS) IF (KCURV.EQ.-1) ED(I)=FF GO TO 300 105 NS=1 C PHASES FOR WEIGHTED F-MAP DO 110 IGP=1,NGP NF=NS+NAG(IGP)-1 IF (NINF(IGP).EQ.4.OR.NINF(IGP).EQ.6.OR.NINF(IGP).EQ.7) GO TO 120 NS=NF+1 110 CONTINUE 120 AC(I)=0.0 BC(I)=0.0 IENT=0 DO 160 J=1,NSYM T=1000.0 DO 150 L=1,3 T=T+FLOAT(I1(L))*TS(L,J) I2(L)=IS(L,1,J)*I1(1) + IS(L,2,J)*I1(2) + IS(L,3,J)*I1(3) 150 CONTINUE CALL SFAC(NS,NF,I2,T,RHO(I),A,B,IENT) IENT=1 AC(I)=AC(I)+A IF (ICENT.EQ.0) BC(I)=BC(I)+B IF (ICENT.EQ.1) AC(I)=AC(I)+A 160 CONTINUE IIG=1 IF (IBGR.EQ.1) IIG=MG+1 FCL=PTS*SQRT(AC(I)*AC(I)+BC(I)*BC(I)) IF (NINF(IGP).EQ.4) EDP(I)=FF/(FCL**2/PTS+FORM*EPS) IF (SC(1).LT.0.0001) SC(1)=1.0 IF (FORM.LT.0.001) ISIM=0 IF (ISIM.EQ.1) 1 XKP(I) = 2.0*SC(1)*EXP(-BT(IIG)*RHO(I))*FCL/FORM/EPS/PTS FC(I)=FCL*EXP(-BT(IIG)*RHO(I)) FO(I)=SC(1)*FO(I) ZC=ZC+FC(I) ZO=ZO+FO(I) ZCG(IG)=ZCG(IG)+FC(I) ZOG(IG)=ZOG(IG)+FO(I) 300 CONTINUE I=600 C CREATE SCRATCH FILE 310 IF (NDIFF.EQ.0) WRITE (8) LH,LK,LL,FO,ID,EW,ED,RHO,FC,EDP C CREATE FILE FOR WEIGHTED FOURIER IF (FO(I).LT.0.0) LH(I)=-1000 IF (NDIFF.EQ.1) THEN WRITE (2) LH,LK,LL,AC,BC,FO,FC,ID,RHO,XKP WRITE (8) LH,LK,LL,AC,BC,FO,FC,ID,RHO,XKP ENDIF RETURN END C ------------------ C BESSEL FUNCTION I1(X)/I0(X) FUNCTION BS10(X) REAL*8 Y,P1,P2,P3,P4,P5,P6,P7,Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9, 1 R1,R2,R3,R4,R5,R6,R7,S1,S2,S3,S4,S5,S6,S7,S8,S9 DATA P1,P2,P3,P4,P5,P6,P7/1.0D0,3.5156229D0,3.0899424D0, 1 1.2067492D0,0.2659732D0,0.360768D-1,0.45813D-2/, 2 Q1,Q2,Q3,Q4,Q5,Q6,Q7,Q8,Q9/0.39894228D0,0.1328592D-1, 3 0.225319D-2,-0.157565D-2,0.916281D-2,-0.2057706D-1, 4 0.2635537D-1,-0.1647633D-1,0.392377D-2/, 5 R1,R2,R3,R4,R5,R6,R7/0.5D0,0.87890594D0,0.51498869D0, 6 0.15084934D0,0.2658733D-1,0.301532D-2,0.32411D-3/, 7 S1,S2,S3,S4,S5,S6,S7,S8,S9/0.39894228D0,-0.3988024D-1, 8 -0.362018D-2,0.163801D-2,-0.1031555D-1,0.2282967D-1, 9 -0.2895312D-1,0.1787654D-1,-0.420059D-2/ IF (ABS(X).LT.3.75) THEN Y=(X/3.75)**2 BS10=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*(R6+Y*R7)))))) 1 /(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))) ELSE Y=3.75/ABS(X) BS10=(S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*(S6+Y*(S7+Y*(S8+Y*S9)))))))) 1 /(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9)))))))) ENDIF RETURN END C----------------------------------------------------------------------- SUBROUTINE INCELL C INPUT UNIT CELL PARAMETERS AND CALCULATE RECIPROCAL PARAMETERS COMMON/UNIT1/ITLE(80),LIST,PI COMMON/XXX/P(6),CX(9),NREF,NB,RHOMAX,MM,EMAX,EN,MZ,ER,ISIM,RHOCUT, 1 RHOMIN,RHOLOW,ENG(8),ERG(8),KNOWN,IPATH,VOLUME CHARACTER ITLE WRITE (6,20) (CX(I),I=1,6) 20 FORMAT(/1X,10HUNIT CELL:,7X,3HA =,F8.3,7X,3HB =,F8.3,7X, 1 3HC =,F8.3/14X,7HALPHA =,F7.2,5X,6HBETA =,F7.2,4X,7HGAMMA =,F7.2) CALL VOL(CX,V) C VOLUME AND RECIPROCAL CELL FUNCTIONS V=1.0/(V*CX(1)*CX(2)*CX(3)) VOLUME=1.0/V P(1)=CX(2)*CX(3)*CX(7)*V P(2)=CX(1)*CX(3)*CX(8)*V P(3)=CX(1)*CX(2)*CX(9)*V P(4)=0.5*P(1)*P(2)*(CX(4)*CX(5)-CX(6))/(CX(7)*CX(8)) P(5)=0.5*P(1)*P(3)*(CX(4)*CX(6)-CX(5))/(CX(7)*CX(9)) P(6)=0.5*P(2)*P(3)*(CX(5)*CX(6)-CX(4))/(CX(8)*CX(9)) DO 40 I=1,3 P(I)=0.25*P(I)*P(I) CX(I+3)=180.0*ATAN2(CX(I+6),CX(I+3))/PI 40 CONTINUE C WRITE (6,50) P C 50 FORMAT(/1X,23H(SIN(THETA)/LAMBDA)**2=, C 1 F8.6,9H * H**2 +,F9.6,9H * K**2 +,F9.6,9H * L**2 +, C 2 /24X,F8.6,10H * H * K +,F8.6,10H * H * L +,F8.6,8H * K * L/) RETURN END ************************************************************************ * * * IIIIIII N N PPPPPP U U TTTTTTT * * I NN N P P U U T * * I N N N P P U U T * * I N N N PPPPPP U U T * * I N N N P U U T * * I N NN P U U T * * IIIIIII N N P UUUUU T * * * * FREE FORMAT DATA INPUT * * VERSION 1998 * ************************************************************************ SUBROUTINE INPUT_P(JUMP,NON,MISS,KMARK,NTOTAL,BTT) COMMON/ATMFACTOR/AL(8),AS(8),BL(8),BS(8),CL(8),CS(8),DL(8),DS(8), 1 EL(8),NW(8),NO(8),NK,NAT,F(9) COMMON/ATMGROUP/NINF(10),NAG(10),X(5000),Y(5000),Z(5000),NZ(5000), 1 Q(5000),U11(5000),U22(5000),U33(5000),U23(5000), 2 U13(5000),U12(5000) COMMON/REFLXIN/LH(600),LK(600),LL(600),FO(600),ID(600),RHO(600), 1 SIG(600),XKP(600) COMMON/RESCALING/ISC,IBGR,MG,MIG(8),SCL,BT(9),SC(9),IP(8,3,5), 1 SCAL(8),RSTT COMMON/SCATFACTOR/GIS(142),GIW(142),NGP,NDIFF,ISOL,X0(3) COMMON/SYMMETRY/IS(3,3,24),TS(3,24),NSYM,PTS,KSYS,ICENT,LATT, 1 IAVE,EXTI COMMON/UNIT1/ ITLE(80),LIST,PI,KCURV,IPAT,IAPA,MOV,MFP,NAU COMMON/XXX/P(6),CX(9),NREF,NB,RHOMAX,MM,EMAX,EN,MZ,ER,ISIM,RHOCUT, 1 RHOMIN,RHOLOW,ENG(8),ERG(8),KNOWN,IPATH,VOLUME DIMENSION POP(10),NA(8),CR(9,10),VECT(8,3),X00(120) DIMENSION ICONV(7500),WTFOM(3),KEYWRD(93) DIMENSION IDIV(20),IB(20),MARKNO(98) CHARACTER N(80),LETT(26),KX(10),KSP,KM,KD,KEQ,KP,KSC,ITERM(4), 1 NGS(26),M(80),ITLE,lett1(26) real cellp(6) DATA LETT/1HA,1HB,1HC,1HD,1HE,1HF,1HG,1HH,1HI,1HJ,1HK,1HL,1HM, 1 1HN,1HO,1HP,1HQ,1HR,1HS,1HT,1HU,1HV,1HW,1HX,1HY,1HZ/ DATA LETT1/1Ha,1Hb,1Hc,1Hd,1He,1Hf,1Hg,1Hh,1Hi,1Hj,1Hk,1Hl,1Hm, 1 1Hn,1Ho,1Hp,1Hq,1Hr,1Hs,1Ht,1Hu,1Hv,1Hw,1Hx,1Hy,1Hz/ DATA KX/1H0,1H1,1H2,1H3,1H4,1H5,1H6,1H7,1H8,1H9/ DATA KSP/1H /,KM/1H-/,KD/1H./ DATA KEQ/1H=/,KP/1H+/,KSC/1H;/ DATA ITERM/1HH,1HK,1HL,1HN/ DATA IOFR/0/,MKREJ/1/,ISTO/0/,NOJOIN/0/,ITAN/0/,NANY/0/,KMIN/60/, 1IPATH/0/,NDET/0/,LSTP/-1/,METAL/0/,PROB/1.01/,ISTAGE/1/, NGEN/0/, 2 NSPEC/0/,NINPUT/0/,IMKK/0/,IFOM/1/,KMAX/5000/,MN/0/,IPUB/0/, 3 ISKIP/0/, CUT1/1.8/,CUT2/1.6/,WTFOM/0.2,1.4,1.4/,NOSET/0/, 6 GRID/0.28/,IWMIN/25/,NRAN/250/,ANGMI/85.0/,ANGMA/145.0/, 7 DMIN/1.1/, DMAX/1.95/,DMUT/2.4/,DM/0.85/,IWT/2/,NAA/0/,MFRN/0/, 8 NPK/0/,NSREQ/0/,DFRG/2.8/,IHAR/0/, NHAR/0/,NHV/0/ DATA KEYWRD 1 /1000, 2000, 3000, 10219, 10321, 11220, 11301, 11309, 11601, 2 12015, 12205, 20321, 20601, 20718, 30512, 30815, 31514, 32220, 3 40520, 40615, 41301, 41309, 41521, 50613, 50914, 51301, 51309, 4 51401, 51619, 51917, 60103, 60119, 60918, 70514, 71809, 71815, 5 80118,110118,110321,111301,111309,111415,120119,120913,120919, 6 130124,130520,130919,131512,131522,132112,140121,140801,140822, 7 141510,141513,141518,141623,141801,141805,142618,151806,151809, 8 151814,160119,160120,160501,160513,160519,160801,161220,161519, 9 161613,161619,161815,161909,180114,180519,180815,181903,181920, * 190301,190520,191109,191512,191605,191620,192015,192320,211423, * 230615,231309,120102/ CALL CCPDPN(7,'SCRA7.TM','SCRATCH','F',80,0) c OPEN(7,FILE='SCRA.TM',FORM='FORMATTED',STATUS='UNKNOWN') C NAME OF THE WHOLE PACKAGE, TO BE PRINTED OUT ONLY WITH JOB REMARKS 10 FORMAT(////// &60H SSSSSSSS AA PPPPPPPPP II / &60H SSSSSSSSSS AAAA PPPPPPPPPP II / &60H SS SS AA AA PP PP II / &60H SS AA AA PP PP II 9999 55555 / &60H SSSSSSSSS AA AA PPPPPPPPPP II 9 9 5 / &60H SSSSSSSSS AA AA PPPPPPPPP II 9 9 5 / &60H SS AAAAAAAAAA PP II === 99999 55555 / &60H SS SS AAAAAAAAAA PP II 9 5 / &60H SSSSSSSSSS AA AA PP II 9 9 5 5 / &60H SSSSSSSS AA AA PP II 9999 5555 / &/////) KUSER3=15000 C SET INITIAL AND DEFAULT VALUES IVAE=0 NHAF=0 MISS=0 DO 12 I=1,120 12 X00(I)=0.0 DO 20 I=1,7500 ICONV(I)=0 IF (I.GT.5000) GO TO 20 X(I)=0.0 Y(I)=0.0 Z(I)=0.0 Q(I)=1.0 U11(I)=0.0 U22(I)=0.0 U33(I)=0.0 U23(I)=0.0 U13(I)=0.0 U12(I)=0.0 IF (I.GT.98) GO TO 20 MARKNO(I)=0 IF (I.GT.20) GO TO 20 IB(I)=0 IDIV(I)=0 IF (I.GT.10) GO TO 20 NINF(I)=0 DO 14 J=1,9 CR(J,I)=0.0 14 CONTINUE IF (I.GT.8) GO TO 20 MIG(I)=1 AL(I)=0.000 AS(I)=0.000 BL(I)=0.000 BS(I)=0.000 CL(I)=0.000 CS(I)=0.000 DL(I)=0.000 DS(I)=0.000 EL(I)=0.000 SCAL(I)=0.0 SC(I)=-1.0 NW(I)=0 NO(I)=0 BT(I)=-200.0 DO 16 J=1,3 VECT(I,J)=0.0 DO 16 K=1,5 IP(I,J,K)=0 16 CONTINUE 20 CONTINUE C DEFAULT UNIT CELL CONTENT: C ATOM --- CARBON. THE NUMBER OF WHICH IS TO BE ESTIMATED C ACCORDING TO THE VOLUME OF UNIT CELL KMARK=0 NK=1 NW(1)=3 NA(1)=-100 JUMP=-1 ITPR=0 NON=0 IND=1 IFAST=0 ISOL=0 MM=0 MG=0 NB=0 KCURV=0 MZ=100 LIST=0 IBGR=0 NTOT=0 NAU=0 ISC=1 NGP=0 NSYM=1 ISIM=1 C -- 3A CUTOFF RHOCUT=0.0278 RHOLOW=0.0 EMAX=5.0 EN=1.2 ER=0.3 SCL=0.5 IPAT=0 MOV=0 MFP=10 IAPA=0 NAPA=0 RSTT=-1.0 NSYM = 1 pts = 1 ICENT = 0 C -- flag latt = 0 cx(1) = 0.0 80 KS=1 MEQ=0 READ (10,90) N 90 FORMAT(80A1) 100 KEY=0 IKW=0 IC=0 IMK=0 IM=1 IGG=0 110 DO 200 I=KS,80 DO 120 K=1,26 IF (LETT(K).EQ.N(I).or.lett1(k).eq.n(i)) GO TO 190 120 CONTINUE IF (IKW.EQ.0) GO TO 200 IF (IKW.LT.3.AND.N(I).EQ.KSP) GO TO 185 IF (IKW.LT.3) GO TO 6000 IF (IC.GT.20) GO TO 200 DO 130 K=1,10 IF (KX(K).EQ.N(I)) GO TO 180 130 CONTINUE IF (N(I).EQ.KD) GO TO 170 IF (N(I).EQ.KM) GO TO 150 IM=1 IF (N(I).EQ.KSP) GO TO 160 IF (N(I).EQ.KSC) GO TO 160 IF (N(I).EQ.KEQ) GO TO 140 IF (N(I).EQ.KP) GO TO 167 GO TO 6000 140 MEQ=MEQ+1 GO TO 167 150 IM=-1 160 IF (IMK.LE.0) IC=IC+1 IF (IC.GT.20) GO TO 200 IF (IC.LE.2) GO TO 167 C (GRO) INDEX GROUP SPECIFICATION IF (KEY.NE.71815) GO TO 167 DO 162 KK=1,4 IF (ITERM(KK).EQ.N(I-1)) GO TO 165 162 CONTINUE IP(MI,MJ,5)=IB(2) 165 IC=2 167 IB(IC)=0 IDIV(IC)=0 IMK=1 GO TO 200 170 IMK=-1 GO TO 200 180 IB(IC)=IM*(10*IABS(IB(IC))+K-1) IF (IMK.EQ.1) IMK=0 IF (IMK.LT.0) IDIV(IC)=IDIV(IC)-1 IF (I.EQ.80) GO TO 181 IF (N(I+1).NE.KSP) GO TO 200 C TEST FOR REMARKS (180513=rem) 181 IF (KEY.EQ.180513) GO TO 300 c (160801=pha, 190520=set, 12015=ato, 131512=mol) IF (KEY.EQ.160801.OR.KEY.EQ.190520.OR.KEY.EQ.12015.OR. 1 KEY.EQ.131512) GO TO 210 c (111415=kno, 161613=ppm, 160513=pem, 160519=pes, 161619=pps c 140822=nhv THREE PARAMETERS (x y z coordinates ) WILL c FOLLOW THOSE KEYWORDS) IF (KEY.EQ.111415.OR.KEY.EQ.161613.OR.KEY.EQ.160513.OR. 1 KEY.EQ.160519.OR.KEY.EQ.161619.OR.KEY.EQ.140822) GO TO 182 NREC=IC IF (IC.NE.20) GO TO 200 IGG=1 KS=I+1 IF (KS.GT.80) THEN KS=1 READ (10,90) N ENDIF IC=0 GO TO 220 182 NREC=IC IF (IC.NE.18) GO TO 200 IGG=1 KS=I+1 IF (KS.GT.80) THEN KS=1 READ (10,90) N ENDIF IC=0 GO TO 220 185 GO TO (6000,6000,2100,2200,2300),IND 190 IF (IC.LE.1.AND.IGG.EQ.1) GO TO 192 IGG=0 IF (IC.GT.0) GO TO 210 GO TO 195 192 KS=I GO TO 100 195 IF (IKW.GE.3) GO TO 200 IKW=IKW+1 KEY=100*KEY+K IF (IKW.LT.3) GO TO 200 C TEST FOR END OF KEYWORD FILE IF (KEY.EQ.51404) GO TO 2800 C (REMARK) IF (KEY.EQ.180513.OR.ITPR.EQ.1) GO TO 197 C PRINT TITLE WRITE (6,360) ITLE ITPR=1 C (SPACE GROUP: SPG) 197 IF (KEY.EQ.191607) GO TO 1630 200 CONTINUE C GRO IF (KEY.EQ.71815) GO TO 80 KS=1 READ (10,90) N GO TO 110 210 KS=I 220 DO 250 II=1,93 IF (KEY.NE.KEYWRD(II)) GO TO 250 GO TO (2150,2250,2350, 400, 420, 430, 440, 450, 465, 2 470, 485, 490, 500, 600, 610, 620, 630, 685, 690, 3 700, 710, 720, 730, 750, 760, 765, 770, 780, 4 790, 800, 810, 820, 830, 840, 850, 860, 887, 5 890, 900, 910, 920, 930, 940, 955, 960, 980, 6 1200,1205,1210,1220,1240,1250,1253,1255,1270, 7 1290,1300,1310,1320,1330,1340,1350,1360,1370, 8 1380,1390,1395,1400,1408,1410,1440,1445,1447, 9 1449,1460,1470,1480,1490,1550,1560,1570,1575, * 1590,1600,1610,1620,1650,1660,1670,1680,1690, * 1700,80), II 250 CONTINUE GO TO 6000 300 WRITE (6,10) ITPR=1 IREM=IB(1) DO 340 I=1,IREM READ (10,90) N WRITE (6,330) N 330 FORMAT(50X,80A1) 340 CONTINUE C PRINT JOB TITLE WRITE (6,360) ITLE 360 FORMAT(1X,'PRELIMINARY PROCESSING OF INPUT ', 1 'DATA'//80A1) GO TO 80 C *** KEYWORD PROCESSING *** C *ABS* ABSFOM WEIGHT 400 WTFOM(1) = FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *ACU* INITIAL CUTOFF OF THE FIRST EARLY FIGURE OF MERIT 420 CUT1 = FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *ALT* LIST A FULL CONVERGENCE MAP 430 LSTP = 1 GO TO 100 C *AMA* MAXIMUM ALLOWED BOND ANGLE 440 ANGMA=FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *AMI* MINIMUM ALLOWED BOND ANGLE 450 ANGMI=FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *APA* AUTOMATIC PATTERSON ANALYSIS 465 NAPA = 1 ISTO=1 GO TO 100 C *ATO* ADDITIONAL ATOMS TO BE CONSIDERED IN "SEARCH" 470 NAA = IB(1) DO 480 I=1,NAA READ (10,90) N WRITE (7,90) N 480 CONTINUE GO TO 80 C *AVE* DATA AVERAGE 485 IAVE=1 GO TO 100 C *BCU* SECOND CUTOFF OF THE EARLY FIGURE OF MERIT 490 CUT2 = FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *BFA* B-FACTOR 500 BTT=0.0 DO 505 J=1,8 IF (J.LT.IC) BT(J)=FLOAT(IB(J))*10.0**IDIV(J) 505 CONTINUE IC1=IC-1 IF (IC1.EQ.1) BTT=BT(1) IF (IC1.GT.1.AND.ISC.EQ.1) ISC=0 JUMP=JUMP+1 GO TO 100 C *BGR* WILSON PLOT FOR EACH REFLECTION GROUPS 600 IBGR=1 NAU=1 GO TO 100 C *CEL* CELL DIMENSIONS 610 IF (IC.NE.7) GO TO 6000 DO 615 J=1,6 IF (NGP.EQ.0) CX(J)=FLOAT(IB(J))*10.0**IDIV(J) IF (NGP.GT.0) CR(J,NGP)=FLOAT(IB(J))*10.0**IDIV(J) 615 CONTINUE IF (NGP.EQ.0) CALL INCELL IF (NA(1).LT.0) NA(1)=VOLUME*0.06023/400. GO TO 100 C *CHO* PHASE-SET NUMBER CHOSEN FOR GENERATING AN E-MAP 620 NOSET = IB(1) GO TO 100 C *CON* CONTENTS OF THE UNIT CELL 630 NK=0 IND=3 GO TO 100 C *CVT* CONVENTIONAL WEIGHTED TANGENT REFINMENT 685 ITAN = 2 GO TO 100 C *DET* NUMBER OF PHASES TO BE DETERMINED 690 NDET = IB(1) GO TO 100 C *DFO* DIFFERENCE FOURIER 700 IND=7 ISTO = 1 GO TO 1710 C *DMA* MAXIMUM ALLOWED BOND LENGTH 710 DMAX=FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *DMI* MINIMUM ALLOWED BOND LENGTH 720 DMIN=FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *DOU* MAXIMUM INTERPEAK DISTANCE TO BE TABULATED 730 DMUT=FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *EFM* EARLY FIGURE OF MERIT IN "PHASE" WILL BE USED 750 IFOM = 0 GO TO 100 C *EIN* E VALUES INSTEAD OF F(OBS) ARE USED AS DIFFRACTION DATA 760 BT(1)=0.0 SCAL(1)=1.0 KCURV=-1 JUMP=1 GO TO 100 C *EMA* MAXIMUM VALUE OF LARGEST E'S PASSED TO "PHASE" 765 EMAX=FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *EMI* MINIMUM VALUE OF LARGEST E'S PASSED TO "PHASE" 770 EN=FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *ENA* ENANTIOMORPH FIXING REFLECTION 780 IEND = IC-1 DO 785 II=1,IEND NINPUT = NINPUT + 1 IF (NINPUT.GT.7500) GO TO 1720 ICONV(NINPUT) = IB(II) * 1000000 + 300 785 CONTINUE GO TO 100 C *EPS* MAXIMUM VALUE OF PSI(ZERO) E'S PASSED TO "PHASE" 790 ER=FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *ESQ* UPPER LIMIT OF FOR 'WEAK' REFLECTION GROUPS 800 SCL=FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *FAC* ANALYTICAL CONSTANTS OF ATOMIC SCATTERING FACTORS 810 IND=5 GO TO 100 C *FAS* RUNNING THE WHOLE JOB IN A FAST WAY 820 IFAST=1 NRAN=100 NAU=1 GO TO 100 C *FIR* RUN THE FIRST PART OF "PHASE" 830 IPATH = 1 GO TO 100 C *GEN* NUMBER OF GENERAL REFLECTIONS IN STARTING SET 840 NGEN = IB(1) GO TO 100 C *GRI* GRID SPACING USED IN "EXFFT" 850 GRID=FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *GRO* INDEX GROUP SPECIFICATION 860 MI=IB(1) MJ=MEQ IF (IB(IC).EQ.0) IB(IC)=IM DO 865 K=1,4 IF (ITERM(K).EQ.N(I)) GO TO 870 865 CONTINUE 870 IF (K.LE.3) MJ=MJ+1 IP(MI,MJ,K)=IB(IC) KS=KS+1 GO TO 110 C *HAR* TEST FOR HARKER ANALYSIS 887 IAPA=IAPA+1 ISTO=1 IF (IB(1).EQ.0) GO TO 100 C IHAR --- NUMBER OF HARKER PEAKS IHAR=IB(1) GO TO 100 C *KAR* KARLE RECYCLING 890 KMARK=IB(1) IND=5 ITAN=1 NRAN=500 IWMIN=45 GO TO 1710 C *KCU* K-CURVE NORMALIZATION 900 KCURV=1 GO TO 100 C *KMA* MAXIMUM KAPPA-VALUE ACCEPTED IN "PHASE" 910 KMAX = INT(100.0*FLOAT(IB(1))*10.0**IDIV(1)) GO TO 100 C *KMI* MINIMUM KAPPA-VALUE ACCEPTED IN "PHASE" 920 KMIN = INT(100.0*FLOAT(IB(1))*10.0**IDIV(1)) GO TO 100 C *KNO* KNOWN PHASES 930 DO 935 II=1,NREC,3 NINPUT = NINPUT + 1 IF (NINPUT.GT.7500) GO TO 1720 IPACK = IB(II) * 1000000 C READ PHASE IPACK = IPACK + IB(II+1) * 1000 C READ WEIGHT ICONV(NINPUT) = IPACK + IB(II+2) 935 CONTINUE IF (IGG.EQ.1) GO TO 110 GO TO 100 C *LAS* RUN THE LAST PART OF PHASE 940 IPATH = 2 NON = 1 GO TO 100 C *LIM* RESOLUTION LIMIT FOR REFLECTIONS TO BE ACCEPTED(ANGSTROM) 955 ANSTR1=FLOAT(IB(1))*10.0**IDIV(1) ANSTR2=FLOAT(IB(2))*10.0**IDIV(2) RHOLOW=1/(4*ANSTR1*ANSTR1) RHOCUT=1/(4*ANSTR2*ANSTR2) GO TO 100 C *LIS* LIST ALL F(OBS) AND E'S 960 LIST=1 GO TO 100 C *MAX* MAXIMUM NUMBER OF PHASE SETS TO BE GENERATED 980 NSREQ = IB(1) GO TO 100 C *MET* NUMBER OF HEAVY ATOMS TO BE CONSIDERED IN "SEARCH" 1200 METAL=IB(1) GO TO 100 C *MIS* COMPENSATION OF THE MISSING REFLECTIONS FOR WILSON STAT. 1205 MISS=1 GO TO 100 C *MOL* SPECIFYING MOLECULAR CONNECTIVITY FOR "SEARCH" 1210 MFRN=IB(1) DO 1215 I=1,MFRN READ (10,90) N WRITE (7,90) N 1215 CONTINUE GO TO 80 C *MOV* REMOVING PATTERSON ORIGIN 1220 MOV=1 GO TO 100 C *MUL* MULTIPLICITY C NUMBER OF TIMES THE RONDOM GROUP APPEARS IN THE UNIT CELL 1240 POP(NGP)=FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *NAU* BY-PASS THE SUBROUTINE "AUTOGP" 1250 NAU=1 GO TO 100 C *NHA* TEST FOR NON-HARKER ANALYSIS 1253 IAPA=IAPA+2 ISTO=1 C NHAR -- NUMBER OF NON-HARKER PEAKS IF (IB(1).EQ.0) GO TO 100 NHAR=IB(1) GO TO 100 C *NHV* NUMBER OF HEAVY ATOMS USED FOR NON-HARKER ANALYSIS 1255 II=NREC/3 IF (II.EQ.0) GO TO 1257 DO 1256 J=1,II NHAF=NHAF+1 IF (NHAF.GT.40) GO TO 6000 K=(NHAF-1)*3 X00(K+1)=FLOAT(IB(K+1))*10.0**IDIV(K+1) X00(K+2)=FLOAT(IB(K+2))*10.0**IDIV(K+2) X00(K+3)=FLOAT(IB(K+3))*10.0**IDIV(K+3) 1256 CONTINUE IF (IGG.EQ.1) GO TO 110 NHV=NHAF GO TO 100 1257 IF (IB(1).GT.40) GO TO 6000 NHV=IB(1) GO TO 100 C *NOJ* NO INTERPRETATION OF PEAK CONNECTIVITY SHOULD BE MADE 1270 NOJOIN=1 GO TO 100 C *NOM* NO MODIFICATION ON RUNNING "PHASE" DISREGARDING WHETHER C THE STRUCTURE CONTAINS PSEUDO-TRANSLATIONAL SYMMETRY 1290 ISTAGE = 0 GO TO 100 C *NOR* NO REJECTING OF 'W-W-W' PHASE RELATIONSHIPS 1300 MKREJ=0 GO TO 100 C *NPW* NUMBER OF POINTS ON WILSON PLOT 1310 NB=IB(1) GO TO 100 C *NRA* NUMBER OF REFLECTIONS IN STARTING SET OF "PHASE" 1320 NRAN = IB(1) GO TO 100 C *NRE* NUMBER OF LARGEST E'S PASSED TO "PHASE" 1330 MM=IB(1) GO TO 100 C *NZR* NUMBER OF PSI-ZERO REFLECTIONS PASSED TO "PHASE" 1340 MZ=IB(1) GO TO 100 C *ORF* FREE OF ORIGIN FIXATION 1350 IOFR=1 GO TO 100 C *ORI* ORIENTED MOLECULAR GROUP 1360 IND=3 GO TO 1710 C *ORN* CODES OF ORIGIN FIXING REFLECTIONS 1370 IEND = IC-1 DO 1375 II=1,IEND NINPUT = NINPUT + 1 IF (NINPUT.GT.7500) GO TO 1720 ICONV(NINPUT) = IB(II) * 1000000 + 200 1375 CONTINUE GO TO 100 C *PAS* BY-PASS "PREPAR" AFTER KEYWORDS HAVE BEEN INPUT 1380 NON=1 GO TO 100 C *PAT* PATTERSON 1390 IPAT=2 IF (N(I).EQ.LETT(6).AND.N(I+1).EQ.LETT(6)) IPAT=1 IF (N(I).EQ.LETT(5).AND.N(I+1).EQ.LETT(5)) IPAT=3 KS=I+2 ISTO=1 GO TO 100 C *PEA* NUMBER OF PEAKS TO BE SEARCHED FOR IN THE DENSITY MAP 1395 NPK=IB(1) GO TO 100 C *PEM* MINIMUM FUNCTION OF PATTERSON-MAP AND E-MAP 1400 MFP = 1 1401 ISTO = 1 II=NREC/3 IF (II.EQ.0.AND.MN.EQ.0) GO TO 1403 IF (II.EQ.0) GO TO 6000 DO 1402 J=1,II MN = MN+1 IF (MN.GT.8) GO TO 6000 DO 1402 JJ=1,3 K=(J-1)*3+JJ VECT(MN,JJ)=FLOAT(IB(K))*10.0**IDIV(K) 1402 CONTINUE IF (IGG.EQ.1) GO TO 110 GO TO 100 1403 MN=IB(1) IF (MN.GT.8) GO TO 6000 CALL CCPDPN(15,'SAPIPKS','UNKNOWN','F',80,0) c OPEN(15,FILE='SAPI95.PKS',FORM='FORMATTED',STATUS='UNKNOWN') READ (15,90) M DO 1405 J=1,MN READ (15,1406) (VECT(J,JJ),JJ=1,3) 1405 CONTINUE 1406 FORMAT(5X,3F10.4) CLOSE (15) GO TO 100 C *PES* SUM FUNCTION OF E-MAP AND PATTERSON MAP 1408 MFP=3 GOTO 1401 C *PHA* STARTING PHASES IN A SINGLE RUN OF TANGENT REFINEMENT 1410 IMKK=IB(1) IMKMAX=KUSER3/2 IF (IMKK .GT. IMKMAX) GO TO 1430 IF (IMKK.EQ.0) GO TO 6000 CALL CCPDPN(9,'KARLE.TM','SCRATCH','F',80,0) c OPEN(UNIT=9,FILE='KARLE.TM',FORM='FORMATTED',STATUS='UNKNOWN') 1415 READ (10,90) N DO 1421 JX=1,80 DO 1420 J=1,26 IF (LETT(J).EQ.N(JX)) GO TO 1425 1420 CONTINUE 1421 CONTINUE WRITE (9,90) N GO TO 1415 1425 KS=1 MEQ=0 c CLOSE (9) GO TO 100 1430 WRITE (6,1435) IMKK,IMKMAX 1435 FORMAT(//44H NUMBER OF STARTING PHASES INPUT TO PHASE =, 1 I5//17X,36HMAXIMUM NUMBER ALLOWED IS KUSER3/4 =,I5) GO TO 6000 C *PLT* PARTIAL LIST OF CONVERGENCE MAP 1440 LSTP = 0 GO TO 100 C *POS* POSITIONED MOLECULAR GROUP 1445 IND=4 GO TO 1710 C *PPM* PATTERSON MINIMUM FUNCTION 1447 MFP = 2 GO TO 1401 C *PPS* PATTERSON SUM FUNCTION 1449 MFP = 4 GO TO 1401 C *PRO* LOWEST PROBABILITY FOR ACCEPTANCE OF SIGMA-1 PHASE 1460 PROB = FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *PSI* WEIGHT OF THE PSIZERO FIGURE OF MERIT 1470 WTFOM(2) = FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *RAN* RANDOMLY POSITIONED MOLECULAR GROUP 1480 IND=2 GO TO 1710 C *RES* WEIGHT OF THE RISDUAL FIGURE OF MERIT 1490 WTFOM(3) = FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *RHO* MAXIMUM VALUE OF (SIN(THETA)/LAMBDA)**2 ACCEPTED 1550 RHOCUT=FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *RSC* RESCALE 1560 IF (ISC.NE.2) ISC=0 WRITE (6,1565) 1565 FORMAT(/1X,24X,30HRESCALING FOR EACH INDEX GROUP) GO TO 100 C *RST* RATIO BETWEEN THE NUMBER OF 'STRONG' REFLECTIONS C AND THAT OF THE TOTAL 1570 RSTT=FLOAT(IB(1))*10.0**IDIV(1) GO TO 100 C *SCA* SCALE FACTORS FOR INDIVIDUAL REFLECTION GROUPS 1575 DO 1577 J=1,8 IF (J.LT.IC) SCAL(J)=FLOAT(IB(J))*10.0**IDIV(J) 1577 CONTINUE IC1=IC-1 IF (IC1.GT.1.AND.ISC.EQ.1) ISC=0 JUMP=JUMP+1 GO TO 100 C *SET* SPECIFYING PHASE SETS TO BE DEVELOPED 1590 IPUB = MIN0(IB(1),KUSER3) IF (IPUB.EQ.0) GO TO 6000 CALL CCPDPN(9,'KARLE.TM','SCRATCH','F',80,0) c OPEN(UNIT=9,FILE='KARLE.TM',FORM='FORMATTED',STATUS='UNKNOWN') 1592 READ (10,90) N DO 1596 JX=1,80 DO 1594 J=1,26 IF (LETT(J).EQ.N(JX)) GO TO 1598 1594 CONTINUE 1596 CONTINUE WRITE (9,90) N GO TO 1592 1598 KS=1 MEQ=0 c CLOSE (9) GO TO 100 C *SKI* SKIP THE FIRST N SETS IN PHASE DEVELOPMENT 1600 ISKIP = IB(1) GO TO 100 C *SOL* SOLVING ENANTIOMORPH AMBIGUITY 1610 IND=8 ISOL=1 NRAN=500 X0(1)=0.0 X0(2)=0.0 X0(3)=0.0 IF (IC.LE.1) GO TO 1710 IF (IC.NE.4) GO TO 6000 X0(1)=FLOAT(IB(1))*10.0**IDIV(1) X0(2)=FLOAT(IB(2))*10.0**IDIV(2) X0(3)=FLOAT(IB(3))*10.0**IDIV(3) GO TO 1710 C *SPE* SPECIAL RESCALING 1620 ISC=2 MG=IB(1) GO TO 100 C *SPG* SET UP LINE FOR SPGR; READ TO SPACE, THEN TO NON-SPACE 1630 KS=MAX0(1,MOD(KS+1,81)) IF (N(KS).NE.KSP) GO TO 1630 1640 KS=MAX0(1,MOD(KS+1,81)) IF (N(KS).EQ.KSP) GO TO 1640 CALL SPGR(N,KS,IERR) IF (IERR.EQ.1) GO TO 6000 GO TO 80 C *SPT* NUMBER OF SPECIALS IN STARTING SET 1650 NSPEC = IB(1) GO TO 100 C *STO* NO RUNNING OF "PHASE" 1660 ISTO = 1 GO TO 100 C *SWT* STATISTICALLY WEIGHTED TANGENT REFINEMENT 1670 ITAN = 1 GO TO 100 C *UNW* UNIT WEIGHT FOR FOURIER OR DEFFERENCE FOURIER 1680 ISIM=0 GO TO 100 C *WFO* WEIGHTED FOURIER 1690 IND=6 ISTO = 1 GO TO 1710 C *WMI* INITIAL WEIGHT FOR THE RANDOM STARTING PHASES 1700 IWMIN = INT(100.0*FLOAT(IB(1))*10.0**IDIV(1)) GO TO 100 1710 NGP=NGP+1 NINF(NGP)=IND IND=4 NAG(NGP)=0 POP(NGP)=1.0 GO TO 100 1720 WRITE (6,1730) 1730 FORMAT (//1X,48H** MORE THAN 250 REFLEXIONS INPUT TO CONVERGE **) GO TO 6000 2100 NK=NK+1 NW(NK)=KEY KEY=1000 IKW=3 GO TO 110 2150 NA(NK)=IB(1) GO TO 100 2200 NAG(NGP)=NAG(NGP)+1 NTOT=NTOT+1 DO 2210 I=1,NK IF (KEY.EQ.NW(I)) GO TO 2230 2210 CONTINUE GO TO 6000 2230 NZ(NTOT)=I KEY=2000 IKW=3 GO TO 110 2250 IF (IC.GT.1) X(NTOT)=FLOAT(IB(1))*10.0**IDIV(1) IF (IC.GT.2) Y(NTOT)=FLOAT(IB(2))*10.0**IDIV(2) IF (IC.GT.3) Z(NTOT)=FLOAT(IB(3))*10.0**IDIV(3) IF (IC.GT.5) Q(NTOT)=FLOAT(IB(5))*10.0**IDIV(5) IF (IC.GT.6) U11(NTOT)=FLOAT(IB(6))*10.0**IDIV(6) IF (IC.GT.7) U22(NTOT)=FLOAT(IB(7))*10.0**IDIV(7) IF (IC.GT.8) U33(NTOT)=FLOAT(IB(8))*10.0**IDIV(8) IF (IC.GT.9) U23(NTOT)=FLOAT(IB(9))*10.0**IDIV(9) IF (IC.GT.10) U13(NTOT)=FLOAT(IB(10))*10.0**IDIV(10) IF (IC.GT.11) U12(NTOT)=FLOAT(IB(11))*10.0**IDIV(11) IF (IC.LE.6) GO TO 100 BTT=0.0 DO 2270 J=1,8 BT(J)=0.0 2270 CONTINUE GO TO 100 2300 DO 2310 JF=1,NK IF (KEY.EQ.NW(JF)) GO TO 2320 2310 CONTINUE GO TO 6000 2320 KEY=3000 IKW=3 GO TO 110 2350 IF (IC.GT.1) AL(JF)=FLOAT(IB(1))*10.0**IDIV(1) IF (IC.GT.2) AS(JF)=FLOAT(IB(2))*10.0**IDIV(2) IF (IC.GT.3) BL(JF)=FLOAT(IB(3))*10.0**IDIV(3) IF (IC.GT.4) BS(JF)=FLOAT(IB(4))*10.0**IDIV(4) IF (IC.GT.5) CL(JF)=FLOAT(IB(5))*10.0**IDIV(5) IF (IC.GT.6) CS(JF)=FLOAT(IB(6))*10.0**IDIV(6) IF (IC.GT.7) DL(JF)=FLOAT(IB(7))*10.0**IDIV(7) IF (IC.GT.8) DS(JF)=FLOAT(IB(8))*10.0**IDIV(8) IF (IC.GT.9) EL(JF)=FLOAT(IB(9))*10.0**IDIV(9) IF (IC.GT.10) GO TO 6000 GO TO 100 C END OF KEYWORD FILE 2800 IF (NON.EQ.1) GO TO 5000 C MARK FOR WILSON STATISTICS WITH KNOWN ATOMIC POSITIONS KNOWN=0 DO 2802 L=1,10 IF (NINF(L).EQ.4) KNOWN=1 2802 CONTINUE C CALCULATE ATOMIC SCATTERING FACTORS CALL ATMCOEF NAT=0 DO 2805 L=1,80 M(L)=KSP 2805 CONTINUE if (na(1).lt.0) then CALL LROPEN (1,'HKLIN',1,IFAIL) IF (IFAIL.EQ.1) CALL CCPERR(1,'Error opening HKLIN') CALL LRCELL(1,CELLP) do i = 1, 6 cx(i) = cellp(i) end do call incell NA(1)=VOLUME*0.036/240. CALL lrclos(1) end if DO 2810 L=1,NK K=NW(L)/100 J=NW(L)-100*K IF (K.GT.0) M(L)=LETT(K) IF (J.GT.0) M(L+40)=LETT(J) NW(L)=NA(L) C NO(L)=INT(AL(L)+BL(L)+CL(L)+DL(L)+EL(L)+0.5) IF (NO(L).GT.1) NAT=NAT+NA(L) MARKNO(NO(L))=L 2810 CONTINUE NRHV=0 NRLT=0 NOLT=0 DO 2820 L=2,98 IF (MARKNO(L).EQ.0) GOTO 2820 IF (FLOAT(NAT-NRLT)/NAT.GT.0.1.OR.FLOAT(L).LT.1.8*NOLT) 1 GOTO 2812 NRHV=NAT-NRLT GOTO 2822 2812 NRLT=NRLT+NA(MARKNO(L)) NOLT=L IF (NAT-NRLT.EQ.0) GOTO 2822 2820 CONTINUE 2822 IF (IHAR.EQ.0.AND.NRHV.GT.0) 1 IHAR=NRHV/PTS-NRHV/(PTS*NSYM*(ICENT+1))+2 IF (IHAR.EQ.0) IHAR=2*(NSYM*(ICENT+1)-1) IF (NHAR.EQ.0.AND.(IAPA.GT.1.AND.IAPA.LT.4)) THEN NHAR=NRHV*(NAT-NRHV)/(PTS*PTS*NSYM*(ICENT+1)) IF (NHAR.EQ.0) NHAR=100 ENDIF WRITE (6,2825) (M(L),M(L+40),NA(L),NO(L),L=1,NK) WRITE (6,2826) (M(L),M(L+40), 1 AL(L),AS(L),BL(L),BS(L),CL(L),CS(L),DL(L),DS(L),EL(L),L=1,NK) 2825 FORMAT(/1X,19HUNIT CELL CONTENTS:/ 1 1X,22X,4HATOM,3X,14HNUMBER IN CELL,2X,13HATOMIC NUMBER/ 2 (1H ,24X,2A1,7X,I5,11X,I3)) 2826 FORMAT(/1X,'SCATTERING FACTOR CONSTANTS:'//7X,'F=', 1 'AA*EXP(-A*RHO)+BB*EXP(-B*RHO)+CC*EXP(-C*RHO)+DD*EXP(-D*RHO)+E'/ 2 /8X,'AA A BB B CC C DD D', 3 ' E'/ (1H ,2A1,9F8.3)) DO 2830 L=1,142 GIS(L)=0.0 2830 CONTINUE IF (NGP.EQ.0) GO TO 2880 NF=0 DO 2870 L=1,NGP NS=NF+1 NF=NF+NAG(L) IF (NINF(L).EQ.2) WRITE (6,2842) 2842 FORMAT(/1X,'RANDOMLY ORIENTED CLUSTER(S) FOR WILSON STATISTICS:') IF (NINF(L).EQ.3) WRITE (6,2843) 2843 FORMAT(/1X,'CORRECTLY ORIENTED CLUSTER(S) FOR WILSON STATISTICS:') IF (NINF(L).EQ.4) WRITE (6,2844) 2844 FORMAT(/1X,'CORRECT POSITIONED CLUSTER(S) FOR WILSON STATISTICS:') IF (NINF(L).EQ.5) WRITE (6,2845) 2845 FORMAT(/1X,'ATOMS FOR KARLE RECYCLING:') IF (NINF(L).EQ.6.AND.ISIM.EQ.1) WRITE (6,2846) IF (NINF(L).EQ.6.AND.ISIM.EQ.0) WRITE (6,2847) 2846 FORMAT(/1X,'ATOMS FOR WEIGHTED FOURIER:') 2847 FORMAT(/1X,'ATOMS FOR UNWEIGHTED FOURIER:') IF (NINF(L).EQ.7.AND.ISIM.EQ.1) WRITE (6,2848) 2848 FORMAT(/1X,'ATOMS FOR WEIGHTED DIFFERENCE FOURIER:') IF (NINF(L).EQ.7.AND.ISIM.EQ.0) WRITE (6,2849) 2849 FORMAT(/1X,'ATOMS FOR UNWEIGHTED DIFFERENCE FOURIER:') IF (NINF(L).EQ.8) WRITE (6,2852) X0(1),X0(2),X0(3) 2852 FORMAT(/1X,'ATOMS FOR SOLVING ENANTIOMORPHOUS AMBIGUITY:'/ 1 3X,'PSEUDO-INVERSE CENTRE:', 1 ' X=',F5.2,' Y=',F5.2,' Z=',F5.2) IF (NGP.GT.1) WRITE (6,2853) L,NAG(L),POP(L) 2853 FORMAT(/3X,'CLUSTER',I4,3X, 1 'NUMBER OF ATOMS =',I3,3X,'MULTIPLICITY =',F5.1) WRITE (6,2854) 2854 FORMAT(/9X,'ATOM',11X,'X',14X,'Y',14X,'Z',9X,'OCCUPANCY') DO 2860 J=NS,NF K=NZ(J) IF (NINF(L).LE.2) NA(K)=NA(K)-INT(POP(L)+0.5) IF (NINF(L).EQ.4.OR.NINF(L).EQ.6.OR.NINF(L).EQ.7) GO TO 2856 WRITE (6,2855) M(K),M(K+40),X(J),Y(J),Z(J) 2855 FORMAT(1H ,9X,2A1,3F15.4) GO TO 2860 2856 NA(K)=NA(K)-INT(NSYM*(ICENT+1)*Q(J)*PTS+0.5) WRITE (6,2858) M(K),M(K+40),X(J),Y(J),Z(J),Q(J) 2858 FORMAT(1H ,9X,2A1,4F15.4) 2860 CONTINUE IF (NGP.EQ.1) WRITE (6,2867) NAG(1) 2867 FORMAT(/1X,'NUMBER OF INPUT ATOMS =',I4) C CALCULATE SPHERICALLY AVERAGED MOLECULAR SCATTERING FACTORS DO 2869 J=1,6 F(J)=CR(J,L) 2869 CONTINUE IF (NINF(L).EQ.2) CALL DEBYE(NS,NF,F,POP(L)) 2870 CONTINUE C CALCULATE WILSON (GIW) AND DEBYE (GIS) SCATTERING FACTORS 2880 DO 2900 L=1,142 T=0.01*FLOAT(L-1) TT=T*T GIW(L)=0.0 DO 2890 J=1,NK FZ=AL(J)*EXP(-AS(J)*TT)+BL(J)*EXP(-BS(J)*TT) 1 +CL(J)*EXP(-CS(J)*TT)+DL(J)*EXP(-DS(J)*TT)+EL(J) GIS(L)=GIS(L)+FZ*FZ*FLOAT(NA(J)) GIW(L)=GIW(L)+FZ*FZ*FLOAT(NW(J)) 2890 CONTINUE 2900 CONTINUE 5000 II=0 IF (NAPA.NE.0) IAPA=4 IF (IAPA.NE.0.AND.IPAT.EQ.0) IPAT=2 IF (MFP.NE.10.AND.IPAT.EQ.0) IPAT=2 IF (MFP.EQ.10.AND.IAPA.EQ.0.AND.IPAT.NE.0.AND.NOJOIN.EQ.0) 1 NOJOIN=1 c CLOSE (10,STATUS='DELETE') c close (10) CALL CCPDPN(34,'PHASKW.TM','SCRATCH','F',80,0) c OPEN(4,FILE='PHASKW.TM',FORM='FORMATTED',STATUS='UNKNOWN') C OUTPUT FOR PHASE WRITE (34,5810) PROB,CUT1,CUT2,WTFOM 5810 FORMAT(5X,3HPRO,7X,3HACU,7X,3HBCU,7X,3HABS,7X,3HPSI,7X,3HRES, 1 /6F10.4) WRITE (34,5820) KMIN,IPATH,ISTAGE,LSTP,NSREQ,NANY,NGEN,NSPEC, 1 IFAST,IFOM,ITAN,ISKIP,IMKK,IPUB,ISTO,IOFR,MKREJ, 2 NDET,KMAX,NRAN,IWMIN,NINPUT 5820 FORMAT(' KMI LAS NOM PLT MAX ANY GEN SPT FAS EFM SWT SKI ITL', 1' SET STO OFR NOR'/17I4/' DET KMA NRA WMI INP'/5I5) IF(NINPUT.GT.0) WRITE (34,5830) (ICONV(I),I=1,NINPUT) 5830 FORMAT(6I12) c CLOSE (4) CALL CCPDPN(20,'EXFFKW.TM','SCRATCH','F',80,0) c OPEN(20,FILE='EXFFKW.TM',FORM='FORMATTED',STATUS='UNKNOWN') C OUTPUT FOR EXFFT MAP AND AAPM IF (IPAT.GT.0.AND.MFP.EQ.10.AND.NPK.EQ.0) 1 NPK=INT(FLOAT(NHAR)*1.1+0.5) WRITE (20,5845) NOSET DO 5842 I=1,NK IF (I.NE.1) GO TO 5840 NAMAX=NO(I) NAMIN=NO(I) GO TO 5842 5840 IF (NO(I).LT.6) GO TO 5842 IF (NO(I).GT.NAMAX) NAMAX=NO(I) IF (NO(I).LT.NAMIN) NAMIN=NO(I) 5842 CONTINUE IF (ISOL.NE.1) GO TO 5844 DO 5843 J=1,3 5843 VECT(1,J)=X0(J) 5844 NALIM=NAMAX*NAMIN/4 WRITE (20,5850) ISOL,MFP,MN,IPAT,NALIM,GRID, 1 ((VECT(I,J),J=1,3),I=1,8) WRITE (20,5838) IAPA,IHAR,NHAR,NHV,NPK,DM,X00 5838 FORMAT(5H APA ,4HHAR ,4HNHAR,4H NHV,4H NPK,8H DM /I4, 1 1X,I3,1X,I3,2X,I2,I4,F8.4,10(/12F10.4)) c CLOSE (20) C OUTPUT FOR SEARCH CALL CCPDPN(11,'SEARCH.TM','SCRATCH','F',80,0) c OPEN(11,FILE='SEARCH.TM',FORM='FORMATTED',STATUS='UNKNOWN') WRITE (11,5860) NAA,MFRN,NPK,METAL,NOJOIN,IAPA ANAT=FLOAT(NAT)/(PTS*FLOAT((ICENT+1)*NSYM)) WRITE (11,5870) ANGMI,ANGMA,DMIN,DMAX,DMUT,DM,DFRG,NK,ANAT DO 5854 L=1,NK WRITE (11,5875) M(L),M(L+40),NW(L),NO(L) 5854 CONTINUE 5845 FORMAT(5HNOSET/I5) 5850 FORMAT(' ISO FMP NVE IPA CUT GRID', 1 /5I4,F8.5,8(/3F10.4)) 5860 FORMAT(38H ATOM MOL PEAK MET NOJ APA /6I6) 5870 FORMAT(53H AMIN AMAX DMIN DMAX DMUT DM DFRG/ 1 2F7.2,5F8.4,I4,F8.2) 5875 FORMAT(2A1,2I6) NTOTAL=4+NK+NAA+MFRN IF (NAA.EQ.0) GO TO 5700 DO 5600 I=1,NAA READ (7,90) N WRITE (11,90) N 5600 CONTINUE 5700 IF (MFRN.EQ.0) GO TO 5900 DO 5800 I=1,MFRN READ (7,90) N WRITE (11,90) N 5800 CONTINUE 5900 CONTINUE c CLOSE (11) CLOSE (7) RETURN 6000 WRITE (6,6010) N 6010 FORMAT(14H ERROR IN LINE,5X,80A1) C CLOSE (6) CLOSE (7) STOP' --- ERROR in PREPARE ---' END C----------------------------------------------------------------------- SUBROUTINE OUTPUT(N,M) C PRINT AND OUTPUT TO FILE FOR REFLEXIONS TO BE USED IN PHASE COMMON/REFLXOUT/IX(30000),EX(30000),FX(30000),SCMK(30000) COMMON/RESCALING/ISC,IBGR,MG,MIG(8),SCL,BT(9),SC(9),IP(8,3,5), 1 SCAL(8) COMMON/UNIT1/ITLE(80),LIST,PI DIMENSION J(3),K(3),L(3),E(3),KODE(3) CHARACTER ITLE IF (N.EQ.M) GO TO 60 WRITE (6,10) 10 FORMAT(1X,3(8H CODE,9H H K L,3X,1HE,2X)) NA=N+1 KK=0 DO 50 JJ=NA,M I=JJ IF (N.NE.0) I=M-JJ+NA KK=KK+1 IG=IABS(MOD(IX(I),32)) IX(I)=IX(I)/32 IND=IX(I) J(KK)=IND/65536 IF (IND.LT.0) J(KK)=J(KK)-1 IND=IND-65536*J(KK) K(KK)=IND/256 L(KK)=IND-256*K(KK)-128 K(KK)=K(KK)-128 SCK=FLOAT(ISIGN(1,MIG(IG)))*SQRT(SCAL(IG)) IF (N.NE.0) GO TO 25 SCMK(I)=SCK WRITE (1,20) J(KK),K(KK),L(KK),EX(I),SCMK(I) 20 FORMAT(3I5,2F10.3) GO TO 28 25 WRITE (1,20) J(KK),K(KK),L(KK),EX(I),SCK 28 E(KK)=EX(I) KODE(KK)=JJ-N IF (KK.NE.3) GO TO 50 WRITE (6,30) (KODE(II),J(II),K(II),L(II),E(II),II=1,3) 30 FORMAT(1X,3(3X,I5,3I3,F6.3)) KK=0 50 CONTINUE IF (KK.EQ.0) GO TO 60 WRITE (6,30) (KODE(II),J(II),K(II),L(II),E(II),II=1,KK) 60 KK=0 D=-1.0 WRITE (1,20) KK,KK,KK,D,D RETURN END C----------------------------------------------------------------------- SUBROUTINE PATT C CALCULATE A AND B FOR PATTERSON FUNCTION COMMON/ATMFACTOR/AL(8),AS(8),BL(8),BS(8),CL(8),CS(8),DL(8),DS(8), 1 EL(8),NW(8),NO(8),NK,NAT,F(9) COMMON/REFLXIN/LH(600),LK(600),LL(600),FO(600),ID(600),RHO(600), 1 SIG(600),XKP(600) COMMON/RESCALING/ISC,IBGR,MG,MIG(8),SCL,BT(9),SC(9),IP(8,3,5), 1 SCAL(8) COMMON/UNIT1/ ITLE(80),LIST,PI,KCURV,IPAT,IAPA,MOV,MFP,NAU COMMON/WILSON/FLGW(30,9),FLGD(30,9),AVR(30,9),DCV(50,9),SLOPE(9), 1 FLGK(9),DEL(9),KS(9),EW(600),ED(600),EDP(600) CHARACTER ITLE REWIND 8 100 READ (8) LH,LK,LL,FO,ID,EW,ED,RHO,SIG DO 200 I=1,600 IF (FO(I).LT.0) GO TO 300 D=EXP(BT(1)*RHO(I)) c D=1.0 IF (IPAT.NE.1) GO TO 80 IF (MOV.EQ.0) EW(I)=FO(I)*FO(I)*D*SC(1) IF (MOV.NE.0) EW(I)=FO(I)*FO(I)*D*SC(1)*(1.0-1/EW(I)) GO TO 150 80 IF (IPAT.EQ.3) GO TO 90 IF (MOV.EQ.0) EW(I)=FO(I)*SQRT(ED(I))*D*SC(1) IF (MOV.NE.0) EW(I)=FO(I)*SQRT(ABS(ED(I)-1))*D*SC(1) GO TO 150 90 IF (MOV.EQ.0) EW(I)=ED(I) IF (MOV.NE.0) EW(I)=ABS(ED(I)-1) 150 ED(I)=0.0 200 CONTINUE WRITE (2) LH,LK,LL,EW,ED,RHO,RHO,RHO,RHO GO TO 100 300 LH(I)=-1000 WRITE (2) LH,LK,LL,EW,ED,RHO,RHO,RHO,RHO RETURN END C----------------------------------------------------------------------- SUBROUTINE RECYC(NS,NF,KMARK) C PHASE CALCULATION FOR WEIGHTED KARLE RECYCLING COMMON/ATMFACTOR/AL(8),AS(8),BL(8),BS(8),CL(8),CS(8),DL(8),DS(8), 1 EL(8),NW(8),NO(8),NK,NAT,F(9) COMMON/ATMGROUP/NINF(10),NAG(10),X(5000),Y(5000),Z(5000),NZ(5000), 1 Q(5000),U11(5000),U22(5000),U33(5000),U23(5000), 2 U13(5000),U12(5000) COMMON/REFLXIN/LH(600),LK(600),LL(600),FO(600),ID(600),RHO(600), 1 SIG(600),XKP(600) COMMON/REFLXOUT/IX(30000),EX(30000),FX(30000),SCMK(30000) COMMON/RESCALING/ISC,IBGR,MG,MIG(8),SCL,BT(9),SC(9),IP(8,3,5), 1 SCAL(8) COMMON/SYMMETRY/IS(3,3,24),TS(3,24),NSYM,PTS,KSYS,ICENT,LATT COMMON/UNIT1/ ITLE(80),LIST,PI,KCURV,IPAT,MOE,MFP,NAU COMMON/WILSON/FLGW(30,9),FLGD(30,9),AVR(30,9),DCV(50,9),SLOPE(9), 1 FLGK(9),DEL(9),KS(9),EW(600),ED(600),EDP(600) COMMON/XXX/P(6),CX(9),NREF,NB,RHOMAX,MM,EMAX,EN,MZ,ER,ISIM,RHOCUT 1 ,RHOMIN,RHOLOW DIMENSION I1(3),I2(3),IFAZ(15000),VECT(8,3) DIMENSION WTFOM(3) CHARACTER ITLE C I1/I0 (BESSEL FUNCTIONS) VEC(U)=U*(U+0.4807)/((U+0.8636)*U+1.3943) CALL CCPDPN(4,'KARLE.TM','SCRATCH','F',80,0) c OPEN(4,FILE='KARLE.TM',FORM='FORMATTED',STATUS='UNKNOWN') C CALCULATE PKARL (RATIO OF KNOWN TO COMPLETE STRUCTURE) KARL=0 MARK=-1 DO 150 I=NS,NF K=NZ(I) KARL=KARL+NO(K)*NO(K) 150 CONTINUE KTOT=0 DO 160 I=1,NK KTOT=KTOT+NW(I)*NO(I)*NO(I) 160 CONTINUE PKARL=FLOAT(KARL*NSYM*(ICENT+1))*PTS/FLOAT(KTOT) IF (PKARL.LT.0.25) MARK=1 WRITE (6,180) ITLE,PKARL 180 FORMAT (//1X,76(1H+)//80A1//20X,'PHASE CALCULATION FOR KARLE ', 1'RECYCLING'//19X,39HRATIO OF KNOWN TO COMPLETE STRUCTURE IS,F6.2) IF (PKARL.LT.0.25) PKARL=0.25 IF (PKARL.GT.0.60) PKARL=0.60 C CALCULATE PHASES AND SELECT REFLEXIONS FOR INPUT TO PHASE REWIND 8 SUMF=0.0 SUMC=0.0 200 READ (8) LH,LK,LL,FO,ID,EW,ED,RHO,SIG DO 300 I=1,600 IF (FO(I).LT.0.0) GO TO 310 IIP=65536*LH(I)+256*(LK(I)+128)+LL(I)+128 IG=MOD(ID(I),100) DO 210 IM=1,MM IF (IX(IM).EQ.IIP) GO TO 220 210 CONTINUE GO TO 300 220 SUMF=SUMF+FX(IM) I1(1)=LH(I) I1(2)=LK(I) I1(3)=LL(I) REL=0.0 RIM=0.0 IENT=0 DO 260 J=1,NSYM T=1000.0 DO 250 L=1,3 T=T+FLOAT(I1(L))*TS(L,J) I2(L)=IS(L,1,J)*I1(1) + IS(L,2,J)*I1(2) + IS(L,3,J)*I1(3) 250 CONTINUE CALL SFAC(NS,NF,I2,T,RHO(I),A,B,IENT) IENT=1 REL=REL+A IF (ICENT.EQ.0) RIM=RIM+B IF (ICENT.EQ.1) REL=REL+A 260 CONTINUE E=EX(IM) EX(IM)=PTS*SQRT((REL*REL+RIM*RIM)*EXP(-BT(IG)*RHO(I))) SUMC=SUMC+EX(IM) IFAZ(IM)=0 C CALCULATE SIM WEIGHT XKAP=2.0*E*E*EX(IM)/((1.0-PKARL)*FX(IM)) IF (XKAP.GT.4.0) GO TO 280 IF (XKAP.LT.2.4.OR.E.LT.1.5.OR.EX(IM).LT.PKARL*FX(IM)) GO TO 300 280 IF (ICENT.EQ.0) EX(IM+15000) = VEC(XKAP) IF (ICENT.EQ.1) EX(IM+15000) = TANH(XKAP/2.0) FAZE=(180.0/PI)*ATAN2(RIM,REL)+360.0 IFAZ(IM)=MOD(INT(FAZE+0.5),360) IF (IFAZ(IM).EQ.0) IFAZ(IM)=360 300 CONTINUE GO TO 200 310 SCALE=SUMC/SUMF RS=0.0 NREC=0 DO 350 I=1,MM RS=RS+ABS(SCALE*FX(I)-EX(I)) IF (IFAZ(I).LE.0) GO TO 350 NREC=NREC+1 IX(NREC)=I IX(NREC+15000)=IFAZ(I) FX(NREC+15000)=EX(I+15000) 350 CONTINUE RS=100.0*RS/SUMC IF (NREC.GT.7500) NREC=7500 WRITE (6,360) SCALE,RS,NREC 360 FORMAT (/19X,7HSCALE =,F7.3,14X,10HR-FACTOR =,F7.2//20X, 1 38HNUMBER OF REFLEXIONS PASSED TO PHASE =,I5//1X,76(1H+)) WRITE (4,370) ITLE,NREC 370 FORMAT(80A1/7HPHASES ,I4) WRITE (4,380) (IX(I),I=1,NREC) 380 FORMAT(15I5) WRITE (4,380) (MARK,I=1,NREC) WRITE (4,390) (FX(I+15000),I=1,NREC) 390 FORMAT(15F5.2) IF (KMARK.EQ.0) MARK=555 IF (KMARK.NE.0) MARK=-1 WRITE (4,380) (IX(I+15000),I=1,NREC),MARK c CLOSE (4) RETURN END C----------------------------------------------------------------------- SUBROUTINE RESCA(KSYS,JUMP) C INDEX GROUP RESCALING COMMON/REFLXIN/LH(600),LK(600),LL(600),FO(600),ID(600),RHO(600), 1 SIG(600),XKP(600) COMMON/RESCALING/ISC,IBGR,NN,MIG(8),SCL,BT(9),SC(9),IP(8,3,5), 1 SCAL(8) COMMON/UNIT1/ITLE(80),LIST,PI,KCURV COMMON/WILSON/FLGW(30,9),FLGD(30,9),AVR(30,9),DCV(50,9),SLOPE(9), 1 FLGK(9),DEL(9),KS(9),EW(600),ED(600),EDP(600) COMMON/XXX/P(6),CX(9),NREF,NB,RHOMAX,MM,EMAX,EN,MZ,ER,ISIM,RHOCUT, 1 RHOMIN,RHOLOW,ENG(8),ERG(8),KNOWN,IPATH,VOLUME DIMENSION NG(8),SCS(8) CHARACTER ITLE IF (IBGR.EQ.1.AND.ISC.EQ.1) ISC=0 IF (JUMP.GE.0) KCURV=0 TOT=0.0 IIG=1 IF (IBGR.EQ.1) IIG=NN+1 NW=0 DO 10 I=1,8 SCS(I)=0.0 NG(I)=0 10 CONTINUE IF (KCURV.EQ.0) WRITE (6,80) IF (KCURV.EQ.1) WRITE (6,90) 80 FORMAT(/1X,12X,'***** NORMALIZATION BY LEAST SQUARES STRAIGHT' 1 ,'LINE *****'/) 90 FORMAT(/1X,12X,'***** NORMALIZATION BY K CURVE OR DEBYE CURVE' 1 ,' *****'/) REWIND 8 100 READ (8) LH,LK,LL,FO,ID,EW,ED,RHO,SIG,EDP DO 250 I=1,600 IF (FO(I).LT.0.0) GO TO 300 IG=MOD(ID(I),100) IF (KCURV.EQ.0) ESQ=ED(I)*EXP(BT(IG)*RHO(I)) IF (KCURV.EQ.1) CALL CURVK(ESQ,RHO(I),ED(I),IIG) C UNPACK SYMMETRY FUNCTIONS MULT=ID(I)/10000 TMUL=FLOAT(MULT) TOT=TOT+ESQ*TMUL NW=NW+MULT SCS(IG)=SCS(IG)+ESQ*TMUL NG(IG)=NG(IG)+MULT 250 CONTINUE GO TO 100 300 TOT=TOT/FLOAT(NW) C NN=8 (PARITY GROUPS) FOR TRICLINIC, MONOCLINIC AND ORTHORHOMBIC C NN=6 (MODIFIED PARITY GROUPS) FOR TETRAGONAL C NN=6 (INDEX GROUPS) IN OTHER SYSTEMS C NN=ANY NUMBER FROM 1 TO 8 FOR SPECIAL RESCALING DO 310 I=1,NN IF (NG(I).GT.0) SCS(I)=SCS(I)/(FLOAT(NG(I))*TOT) 310 CONTINUE C FIND OUT AND GIVE NEGATIVE MARK TO C INDEX GROUP HAVING INTENSITIES SYSTEMATICALLY WEAK 315 DO 320 I=1,NN IF (SCS(I).GT.0.00001.AND.SCS(I).LT.0.7) MIG(I)=2 IF (SCS(I).GT.0.00001.AND.SCS(I).LT.SCL) MIG(I)=-1 320 CONTINUE 325 WRITE (6,330) 330 FORMAT(/1X,6X,'AVERAGE E**2 ACCORDING TO APPROPRIATE INDEX ', 1 'GROUP BEFORE RESCALING') IF (ISC.EQ.2) GO TO 382 IF (KSYS.LE.3) WRITE (6,335) 335 FORMAT(/1X,33X,13HPARITY GROUPS/1X,13X,3HALL,4X,3HEEE,4X,3HOEE, 1 4X,3HEOE,4X,3HOOE,4X,3HEEO,4X,3HOEO,4X,3HEOO,4X,3HOOO) IF (KSYS.EQ.4) WRITE (6,340) 340 FORMAT(/1X,28X,22HMODIFIED PARITY GROUPS/1X,13X,3HALL,4X,3HEEE, 1 8X,7HEOE,OEE,8X,3HOOE,4X,3HEEO,8X,7HEOO,OEO,7X,3HOOO) IF (KSYS.LE.4) GO TO 365 WRITE (6,345) 345 FORMAT(/1X,25HINDEX GROUPS DIVIDED ON -) IF (KSYS.LE.6) WRITE (6,350) 350 FORMAT(10X,11H1) MOD(H,3),4X,11H2) MOD(K,3),4X, 1 13H3) MOD(H+K,3),4X,11H4) MOD(L,2)) IF (KSYS.GE.7) WRITE (6,355) 355 FORMAT(10X,13H1) MOD(H-L,3),4X,13H2) MOD(K-L,3),4X, 1 13H3) MOD(H-K,3),4X,15H4) MOD(H+K+L,2)) WRITE (6,360) 360 FORMAT(/1X,15X,18HE - ZERO REMAINDER,4X,22HO - NON-ZERO REMAINDER/ 1 /1X,14X,3HALL,3X,4HOOOE,1X,9HOOEE,OEOE,1X,4HEEEE,1X,4HOOOO, 2 1X,9HOOEO,OEOO,1X,4HEEEO/1H ,27X,4HEOOE,17X,4HEOOO) 365 TT=1.0000 IF(KSYS.EQ.4) THEN WRITE(6,366) TT,(SCS(I),I=1,NN) WRITE(6,368) NW,(NG(I),I=1,NN) ELSE WRITE(6,370) TT,(SCS(I),I=1,NN) WRITE(6,380) NW,(NG(I),I=1,NN) ENDIF 366 FORMAT(1H ,4X,4HE**2,2X,F7.3,F8.3,2F10.3,F8.3,2F10.3) 368 FORMAT(1H ,3X,6HNUMBER,1X,I7,I8,2I10,I8,2I10) 370 FORMAT(1H ,4X,4HE**2,2X,9F7.3) 380 FORMAT(1H ,3X,6HNUMBER,1X,9I7) GO TO 390 382 WRITE (6,383) 383 FORMAT(/1X,29X,20H INDEX GROUP E**2/) WRITE (6,385) (I,SCS(I),I=1,NN) 385 FORMAT(30X,I8,F12.3) 390 IF (JUMP.EQ.1) GO TO 500 TOT=1.0/TOT C TO DETERMINE IF INDEX GROUP RESCALING IS NEEDED. IF (ISC.NE.1) GO TO 398 DO 395 I=1,NN IF (MIG(I).EQ.-1.OR.MIG(I).EQ.2) ISC=0 IF (MIG(I).EQ.-1.OR.MIG(I).EQ.2) GO TO 398 395 CONTINUE WRITE (6,397) 397 FORMAT(/9X,37H*** NO RESCALING FOR INDIVIDUAL INDEX, 1 23H GROUPS IS REQUIRED ***,/) 398 DO 400 I=1,NN IF (NG(I).EQ.0) GO TO 399 SCAL(I)=TOT/SCS(I) 399 IF (ISC.EQ.1) SCAL(I)=TOT 400 CONTINUE 500 IF (ISC.EQ.1) GO TO 541 WRITE (6,540) SCL 540 FORMAT(/1X,10X,44H*** RESCALING FOR INDIVIDUAL INDEX GROUPS IS, 1 13H REQUIRED ***,/9X,34HREFLECTIONS IN INDEX GROUPS HAVING, 2 23H AVERAGE E**2 LESS THAN,F6.3/9X,23HARE TO BE RECOGNIZED AS, 3 34H 'SYSTEMATICALLY WEAK REFLECTIONS') 541 IF (JUMP.EQ.-1) WRITE (6,542) 542 FORMAT(/1X,11X,46HTEMPERATURE AND SCALING FACTORS DERIVED BY THE, 1 8H PROGRAM) IF (JUMP.EQ.0) WRITE (6,543) 543 FORMAT(/1X,20X,37HSCALING FACTOR DERIVED BY THE PROGRAM/ 1 1H ,12X,53H WITH THE TEMPERATURE FACTOR(BT) SUPPLIED BY THE USER) IF (JUMP.EQ.1) WRITE (6,544) 544 FORMAT(/1X,6X,41HTEMPERATURE FACTOR(BT) AND SCALING FACTOR, 1 21H SUPPLIED BY THE USER) WRITE (6,545) 545 FORMAT(1H ,16X,46H--- EXP{(-2*BT)*RHO}*FCAL**2=SCALE*FOBS**2 ---/ 1 1H ,24X,5HGROUP,5X,4H2*BT,6X,5HSCALE) IF (ISC.NE.1) WRITE (6,546) (I,BT(I),SCAL(I),I=1,NN) 546 FORMAT(1H ,27X,I1,2F11.4) IF (ISC.EQ.1) WRITE (6,548) BT(1),SCAL(1) 548 FORMAT(1H ,25X,3HALL,F10.4,F11.4) RETURN END C----------------------------------------------------------------------- SUBROUTINE RFAC(ICENT,NINF) C R-FACTOR CALCULATION COMMON/REFLXIN/LH(600),LK(600),LL(600),FO(600),ID(600),RHO(600), 1 FC(600),XKP(600) COMMON/RESCALING/ISC,IBGR,MG COMMON/XXX/P(6),CX(9),NREF,NB,RHOMAX,MM,EMAX,EN,MZ,ER,ISIM,RHOCUT, 1 RHOMIN,RHOLOW,SUMFCG(8),SUMFOG(8),KNOWN,IPATH,VOLUME,SUMFC,SUMFO DIMENSION RG(8),AO(600),BO(600),SCALEG(8),NINF(10) CALL CCPDPN(26,'FCOUT','UNKNOWN','F',80,0) c OPEN(26,FILE='FC.OUT',FORM='FORMATTED',STATUS='UNKNOWN') WRITE(26,20) 20 FORMAT(' H K L Fc PHAc') RTOD=45.0/ATAN(1.0) REWIND 2 READ (2) READ (2) R=0.0 SCALE=SUMFC/SUMFO DO 50 I=1,8 RG(I)=0.0 SCALEG(I)=0.0 IF(SUMFOG(I).GT.0.01) SCALEG(I)=SUMFCG(I)/SUMFOG(I) 50 CONTINUE REWIND 8 100 READ (8) LH,LK,LL,AO,BO,FO,FC,ID,RHO,XKP DO 200 I=1,600 IF (FO(I).LT.0.0) GO TO 250 IG=MOD(ID(I),100) RG(IG)=RG(IG)+ABS(FC(I)-FO(I)*SCALEG(IG)) FO(I)=SCALE*FO(I) R=R+ABS(FC(I)-FO(I)) C SIM WEIGHTING SCHEME FOR FOURIER WATE=1.0 IF (ISIM.EQ.0) GO TO 110 IF (ICENT.EQ.0) WATE=BS10(FO(I)*XKP(I)) IF (ICENT.EQ.1) WATE=TANH(FO(I)*XKP(I)/2.0) 110 COEFF=0.0 IF (FC(I).LT.0.00001) GO TO 120 IF (NINF(1).EQ.6) COEFF=FO(I)*WATE/FC(I) IF (NINF(1).EQ.7) COEFF=(FO(I)*WATE-FC(I))/FC(I) 120 AO(I)=AO(I)*COEFF BO(I)=BO(I)*COEFF PHAC=RTOD*ATAN2(BO(I),AO(I))+360.0 PHAC=AMOD(PHAC,360.0) IF(ABS(PHAC-360.0).LE.0.1) PHAC=0.0 FCS=FC(I)/SCALE WRITE(26,150) LH(I),LK(I),LL(I),FCS,PHAC 150 FORMAT(I5,2I4,F12.2,F10.2) 200 CONTINUE WRITE (2) LH,LK,LL,AO,BO,FO,FC,ID,RHO GO TO 100 250 R=100.0*R/SUMFC WRITE (2) LH,LK,LL,AO,BO,FO,FC,ID,RHO WRITE (6,280) R,NREF 280 FORMAT(/1X,14X,19HOVER ALL R-FACTOR =,F7.2,7H % FOR,I6,2X, 1 10HREFLEXIONS) DO 290 I=1,MG IF (SUMFCG(I).LT.0.001) GO TO 286 RG(I)=100.0*RG(I)/SUMFCG(I) GO TO 290 286 RG(I)=0.0 290 CONTINUE WRITE (6,300) 300 FORMAT(/1X,20X,37HR-FACTORS FOR INDIVIDUAL INDEX GROUPS/ 1 1H ,30X,16HGROUP R-FACTOR/) WRITE (6,310) (I,RG(I),I=1,MG) 310 FORMAT(1H ,27X,I5,4X,F7.2,' %') CLOSE (26) RETURN END C----------------------------------------------------------------------- SUBROUTINE SFAC(NS,NF,L,T,RHO,A,B,IENT) C STRUCTURE FACTOR CALCULATION COMMON/ATMFACTOR/AL(8),AS(8),BL(8),BS(8),CL(8),CS(8),DL(8),DS(8), 1 EL(8),NW(8),NO(8),NK,NAT,F(9) COMMON/ATMGROUP/NINF(10),NAG(10),X(5000),Y(5000),Z(5000),NZ(5000), 1 Q(5000),U11(5000),U22(5000),U33(5000),U23(5000), 2 U13(5000),U12(5000) COMMON/SINETABLE/SINT(450) COMMON/XXX/P(6),CX(9),NREF,NB,RHOMAX,MM,EMAX,EN,MZ,ER,ISIM,RHOCUT 1 ,RHOMIN,RHOLOW DIMENSION L(3),COST(360) EQUIVALENCE (SINT(91),COST(1)) A=0.0 B=0.0 IF (IENT.EQ.1) GO TO 30 DO 10 I=1,NK F(I)=AL(I)*EXP(-AS(I)*RHO)+BL(I)*EXP(-BS(I)*RHO)+CL(I)*EXP(-CS(I) 1 *RHO)+DL(I)*EXP(-DS(I)*RHO)+EL(I) 10 CONTINUE 30 HJ=FLOAT(L(1)) HK=FLOAT(L(2)) HL=FLOAT(L(3)) DO 50 I=NS,NF N=NZ(I) ARG=AMOD(HJ*X(I)+HK*Y(I)+HL*Z(I)+T,1.0) IARG=INT(360.0*ARG+0.5)+1 IF (IARG.EQ.361) IARG=1 IF (U11(I).GT.0.00001) GO TO 40 FJ=F(N) GO TO 45 40 UU = U11(I)*HJ*HJ*P(1)+U22(I)*HK*HK*P(2)+U33(I)*HL*HL*P(3) * +U12(I)*HJ*HK*P(4)+U13(I)*HJ*HL*P(5)+U23(I)*HK*HL*P(6) CT98 IF (U22(I).LT.0.00001) UU = U11(I)*RHO CT98 FJ=F(N)*EXP(-8*3.1416*3.1416*UU) IF (U22(I).GE.0.00001) UU = 8.0*3.1416*3.1416*UU IF (U22(I).LT.0.00001) UU = U11(I)*RHO FJ=F(N)*EXP(-UU) 45 A=A+FJ*Q(I)*COST(IARG) B=B+FJ*Q(I)*SINT(IARG) 50 CONTINUE RETURN END C----------------------------------------------------------------------- SUBROUTINE SOLENA(NS,NF,MMS) C PHASE CALCULATION FOR SOLVING ENANTIOMORPHOUS AMBIGUITY COMMON/ATMFACTOR/AL(8),AS(8),BL(8),BS(8),CL(8),CS(8),DL(8),DS(8), 1 EL(8),NW(8),NO(8),NK,NAT,F(9) COMMON/ATMGROUP/NINF(10),NAG(10),X(5000),Y(5000),Z(5000),NZ(5000), 1 Q(5000),U11(5000),U22(5000),U33(5000),U23(5000), 2 U13(5000),U12(5000) COMMON/REFLXIN/LH(600),LK(600),LL(600),FO(600),ID(600),RHO(600), 1 SIG(600),XKP(600) COMMON/REFLXOUT/IX(30000),EX(30000),FX(30000) COMMON/RESCALING/ISC,IBGR,MG,MIG(8),SCL,BT(9),SC(9),IP(8,3,5), 1 SCAL(8) COMMON/SCATFACTOR/GIS(142),GIW(142),NGP,NDIFF,ISOL,X0(3) COMMON/SYMMETRY/IS(3,3,24),TS(3,24),NSYM,PTS,KSYS,ICENT,LATT COMMON/UNIT1/ ITLE(80),LIST,PI,KCURV,IPAT,MOE,MFP,NAU COMMON/WILSON/FLGW(30,9),FLGD(30,9),AVR(30,9),DCV(50,9),SLOPE(9), 1 FLGK(9),DEL(9),KS(9),EW(600),ED(600),EDP(600) COMMON/XXX/P(6),CX(9),NREF,NB,RHOMAX,MM,EMAX,EN,MZ,ER,ISIM,RHOCUT 1 ,RHOMIN,RHOLOW DIMENSION I1(3),I2(3),FC(15000),IMARK(15000),WT(15000), 1 IFAZ(15000),IDEL(15000) CHARACTER ITLE C I1/I0 (BESSEL FUNCTIONS) VEC(U)=U*(U+0.4807)/((U+0.8636)*U+1.3943) CALL CCPDPN(4,'KARLE.TM','SCRATCH','F',80,0) c OPEN(4,FILE='KARLE.TM',FORM='FORMATTED',STATUS='UNKNOWN') SCALEF=SQRT(SC(1)) C CALCULATE PKSOL (RATIO OF KNOWN TO COMPLETE STRUCTURE) KSOL=0 KSW=0 IIG=1 IF (IBGR.EQ.1) IIG=MG+1 IF (NF-NS.GT.5) KSW=1 DO 150 I=NS,NF K=NZ(I) KSOL=KSOL+NO(K)*NO(K) 150 CONTINUE KTOT=0 DO 160 I=1,NK KTOT=KTOT+NW(I)*NO(I)*NO(I) 160 CONTINUE PKSOL=FLOAT(KSOL*NSYM*(ICENT+1))*PTS/FLOAT(KTOT) WRITE (6,180) ITLE,PKSOL 180 FORMAT (//1X,76(1H+)//80A1//24X,'SOLVING ENANTIOMORPHOUS AMBIG', 1'UITY'//18X,39HRATIO OF KNOWN TO COMPLETE STRUCTURE IS,F6.2) C CALCULATE PHASES AND DELTA PHASES BY SIM FORMULA REWIND 8 SUMF=0.0 SUMC=0.0 200 READ (8) LH,LK,LL,FO,ID,EW,ED,RHO,SIG DO 300 I=1,600 IF (FO(I).LT.0.0) GO TO 310 IIP=65536*LH(I)+256*(LK(I)+128)+LL(I)+128 DO 210 IM=1,MM IF (IX(IM)/32.EQ.IIP) GO TO 220 210 CONTINUE GO TO 300 220 SUMF=SUMF+FX(IM) I1(1)=LH(I) I1(2)=LK(I) I1(3)=LL(I) REL=0.0 RIM=0.0 IENT=0 DO 260 J=1,NSYM T=1000.0 DO 250 L=1,3 T=T+FLOAT(I1(L))*(TS(L,J)-X0(L)) I2(L)=IS(L,1,J)*I1(1) + IS(L,2,J)*I1(2) + IS(L,3,J)*I1(3) 250 CONTINUE CALL SFAC(NS,NF,I2,T,RHO(I),A,B,IENT) IENT=1 REL=REL+A RIM=RIM+B 260 CONTINUE E=EX(IM) FC(IM)=PTS*SQRT((REL*REL+RIM*RIM)*EXP(-BT(IIG)*RHO(I))) 265 SUMC=SUMC+FC(IM) IF (KSW.EQ.0) GO TO 280 C CALCULATE DELTA PHASE WITH B=SQRT(F*F-A*A) IFAZ(IM)=INT(360.0*AMOD(I1(1)*X0(1)+I1(2)*X0(2)+I1(3)*X0(3)+ 1 100.0,1.0)) FOBS=PKSOL*SCALEF*FO(I)*EXP(0.5*BT(IIG)*RHO(I))/PTS RIM1=FOBS*FOBS-REL*REL IF (RIM1.LT.0.0) RIM1=0.0 RIM1=SQRT(RIM1) DELTA=(180.0/PI)*ATAN2(RIM1,ABS(REL)) WT(IM)=DELTA IF (REL.LT.0.0) IFAZ(IM)=MOD(IFAZ(IM)+180,360) IF (IFAZ(IM).EQ.0) IFAZ(IM)=360 GO TO 300 C CALCULATE SIM WEIGHT 280 XKAP=2.0*E*E*FC(IM)/((1.0-PKSOL)*FX(IM)) FAZE=0.0 IF (FC(IM).LT.0.000001) GO TO 285 FAZE=(180.0/PI)*ATAN2(RIM,REL)+360.0 285 WT(IM)=(180.0/PI)*ABS(ACOS(VEC(XKAP))) IFAZ(IM)=MOD(INT(FAZE+0.5),360) IF (IFAZ(IM).EQ.0) IFAZ(IM)=360 300 CONTINUE GO TO 200 310 SCALE=SUMC/SUMF C CHOOSE MMS REFLECTIONS WITH LARGEST DELTA PHI FROM MM REFLECTIONS IF (MMS.EQ.MM) GO TO 950 NWEAK=0 NBB=0 C FIND THE NUMBERS FOR STRONG AND WEAK REFLECTION DO 500 I=1,MM IMARK(I)=1 IG=IABS(MOD(IX(I),32)) IF (MIG(IG).NE.-1) GO TO 400 NWEAK=NWEAK+1 WT(I)=0.0 IMARK(I)=2 IFAZ(I)=0 GO TO 500 400 IF (WT(I).LT.15.0) GO TO 500 NBB=NBB+1 IMARK(I)=3 500 CONTINUE NSTRON=MM-NWEAK NAA=NSTRON-NBB IF (NSTRON.LT.MMS/2.OR.NWEAK.LT.MMS/2) GO TO 550 NSTRON=MMS/2 NWEAK=MMS-NSTRON GO TO 600 550 IF (NWEAK.LT.MMS/2) NSTRON=MMS-NWEAK IF (NSTRON.LT.MMS/2) NWEAK=MMS-NSTRON 600 IF (NBB.LT.NSTRON*2/3.OR.NAA.LT.NSTRON/3) GO TO 650 NAA=NSTRON/3 NBB=NSTRON-NAA GO TO 700 650 IF (NBB.LT.NSTRON*2/3) NAA=NSTRON-NBB IF (NAA.LT.NSTRON/3) NBB=NSTRON-NAA C DELETE UNSUITABLE REFLECTIONS 700 CALL SORT(WT,EX,IMARK,MM) J=NBB+1 DO 750 I=J,MM IF (WT(I).LT.15.0) WT(I)=0.0 IF (WT(I).GE.15.0) IMARK(I)=0 750 CONTINUE CALL SORT(EX,WT,IMARK,MM) NCAA=0 NCWEAK=0 DO 800 I=1,MM KKSW=IMARK(I)+1 GO TO (800,760,770,800),KKSW 760 IF (NCAA.GT.NAA) GO TO 790 NCAA=NCAA+1 GO TO 800 770 IF (NCWEAK.GT.NWEAK) GO TO 790 NCWEAK=NCWEAK+1 GO TO 800 790 IMARK(I)=0 800 CONTINUE DO 900 I=1,MMS 830 IF (IMARK(I).NE.0) GO TO 880 DO 850 J=I,MM IMARK(J)=IMARK(J+1) IX(J)=IX(J+1) EX(J)=EX(J+1) IFAZ(J)=IFAZ(J+1) WT(J)=WT(J+1) FC(J)=FC(J+1) FX(J)=FX(J+1) 850 CONTINUE GO TO 830 880 IDEL(I)=INT(WT(I)) WT(I)=1.0 900 CONTINUE MM=MMS GO TO 1100 950 DO 1000 I=1,MM IDEL(I)=INT(WT(I)) WT(I)=1.0 1000 CONTINUE 1100 RS=0.0 WRITE (4,1320) 1320 FORMAT(1X,4HCODE,4X,1HH,4X,1HK,4X,1HL,5X,2HFO,5X,2HFC,6X,1HE, & 5X,2HWT,1X,4HPHIK,1X,4HDPHA/) DO 1350 I=1,MM RS=RS+ABS(SCALE*FX(I)-FC(I)) IND=IX(I)/32 IHH=IND/65536 IF (IND.LT.0) IHH=IHH-1 IND=IND-65536*IHH IKK=IND/256 ILL=IND-256*IKK-128 IKK=IKK-128 FOO=FX(I) FCC=FC(I) EO=EX(I) IPHIQ=IFAZ(I) IPHID=IDEL(I) WEIT=WT(I) WRITE (4,1325) I,IHH,IKK,ILL,FOO,FCC,EO,WEIT,IPHIQ,IPHID 1325 FORMAT(4I5,2F7.2,2F7.4,