C C C This code is distributed under the terms and conditions of the C CCP4 licence agreement as `Part 2' (Annex 2) software. C A copy of the CCP4 licence can be obtained by writing to the C CCP4 Secretary, Daresbury Laboratory, Warrington WA4 4AD, UK. C C SUBROUTINE HKON C C---Subroutine for refinement with diagonal approximation INCLUDE 'atom_com.fh' INCLUDE 'pls_incl.fh' INCLUDE 'celsym.fh' INCLUDE 'cif_incl.fh' INCLUDE 'refi_flags.fh' INCLUDE 'agreem.fh' INCLUDE 'const.fh' INCLUDE 'weights.fh' INCLUDE 'vitals.fh' INCLUDE 'tls.fh' C C----These arrays should be controlled dynamicly as they depend C----on coordinates REAL AM,V,DV COMMON /LSQMTX/ AM( QQM),V( QQV),DV( QQV) INTEGER NE,NZZ,NWTB,NSYM_DIST(4,QQD) COMMON /DISTNS/ NE( QQD),NZZ( QQD),NWTB( QQD),NSYM_DIST C C---Local variables INTEGER I_INTER,I_REF_FLAG CHARACTER LINE*80 C C---report about type of refinement and so on C---------------------------------------------------- C------------------------------------------------- C---Initialize IF((NTLS_CYCLE.GT.0.OR.TLS_FILE_FLAG).AND.REFID.NE.'IDEA') THEN UOVER_ATOM = 0.0 CALL TLS_REFINE ENDIF c CALL SET_HKON_FLAGS C c---Save vital parameters CALL SAVE_VITAL C C----Set refinement style IF(REFID.EQ.'REST') THEN GEOM_FLAG = .TRUE. X_RAY_FLAG = .TRUE. ELSEIF(REFID.EQ.'UNRE') THEN GEOM_FLAG = .FALSE. X_RAY_FLAG = .TRUE. ELSEIF(REFID.EQ.'IDEA') THEN GEOM_FLAG = .TRUE. X_RAY_FLAG = .FALSE. ELSE WRITE(LINE,'(A,A)') 'There is no such option ',REFID CALL ERRWRT(-1,LINE) CALL ERRWRT(1,' False option') ENDIF IF(.NOT.X_RAY_FLAG) THEN C C---If we are not going to refine ITEMP = 0 NOCC = 0 ENDIF IF(.NOT.GEOM_FLAG) THEN NDIS = 0 NBIS = 0 NVDW = 0 ENDIF C c---Add these variables to vitals.fh NVPOS = 0 NMPOS = 0 NVOCC = 0 NMOCC = 0 NVTMP = 0 NMTMP = 0 C C---Number of refined parameters. It is useful for SU calculation NREF_PARS = 0 C C---Find number of variables and second derivative matrix elements C---for different type of parameters IF(OCCUP_REF_FLAG)CALL AUTO_DEFINE_OCCUPANCY_GROUPS DO IA=1,N_ATOM I_INTER = ATOM_REF_FLAG(IA)/10 I_REF_FLAG = ATOM_REF_FLAG(IA) - I_INTER*10 IF(ATOM_REF_FLAG(IA).GT.0) THEN NVPOS = NVPOS + 3 NMPOS = NMPOS + 6 IF(I_REF_FLAG.GT.2) NREF_PARS = NREF_PARS + 3 IF(IOCCUP.GT.0) THEN NVOCC = NVOCC + 1 IF(I_REF_FLAG.GT.2) NREF_PARS = NREF_PARS + 1 NMOCC = NMOCC + 1 ENDIF IF(ITEMP.NE.0) THEN IF(U_ANISO(2,IA).EQ.0.0) THEN NVTMP = NVTMP + 1 IF(I_REF_FLAG.GT.2) NREF_PARS = NREF_PARS + 1 NMTMP = NMTMP + 1 ELSE NVTMP = NVTMP + 6 IF(I_REF_FLAG.GT.2) NREF_PARS = NREF_PARS + 6 NMTMP = NMTMP + 21 ENDIF ENDIF ENDIF ENDDO NVALL = NVPOS + NVOCC + NVTMP c C---Begin refinement cycles C C---Memory allocation should be done here IF(SOLVMIN.EQ.'CGMA') THEN CALL CGSPARCE ELSEIF(SOLVMIN.EQ.'CGRA') THEN CONTINUE CALL CGRAD ELSEIF(SOLVMIN(1:4).EQ.'CDIR') THEN CONTINUE CALL CDIR ELSE LINE = ' Such minimization procedure is not ready' CALL ERRWRT(-1,LINE) CALL ERRWRT(1,' Minimization method') ENDIF C cd CALL UGTENV('XYZOUT',NAMEC) cd PATHC = ' ' cd EXTC = ' ' cd MDOC = 0 cd CALL WRITE_ATOMS(MDOC,PATHC,NAMEC,EXTC,IERROR) cd IF(IERROR.GT.0) THEN cd CALL ERRWRT(-1,'Write atoms failed in HKON') cd CALL ERRWRT(1, cd + 'FATAL ERROR. Check if you have specified filename') cd ENDIF cd CALL ADD_SCALE_CONTR CALL RESTORE_VIT RETURN END C C---Following subroutines are methods of minimizations C SUBROUTINE CGSPARCE cd implicit none c c---Minimization of functional by calculation C---sparce matrix and solution of sparce matrix C---linear equation by conjugate gradient method C---See for example: Numerical Recipes in Fortran by Press et al INCLUDE 'atom_com.fh' INCLUDE 'pls_incl.fh' INCLUDE 'cif_incl.fh' INCLUDE 'refi_flags.fh' INCLUDE 'agreem.fh' INCLUDE 'weights.fh' INCLUDE 'vitals.fh' INCLUDE 'monitor.fh' INCLUDE 'occupancy_params.fh' INCLUDE 'save_all_params.fh' INCLUDE 'tls.fh' INTEGER MCOLS PARAMETER (MCOLS = 200) REAL AM,V,DV COMMON /LSQMTX/ AM( QQM),V( QQV),DV( QQV) INTEGER NE,NZZ,NWTB,NSYM_DIST(4,QQD) COMMON /DISTNS/ NE( QQD),NZZ( QQD),NWTB( QQD),NSYM_DIST C REAL VAL,GAMMA COMMON /R_SCRATCH/ VAL(QQDEN) C INTEGER LOOKUP(MCOLS) INTEGER NLPRGI,NLPRGO COMMON /MTZRD/NLPRGI,NLPRGO,LOOKUP EXTERNAL SCALE_SHIFTS_ALL,SCALE_UVAL,MATMUL_UVAL, + SCALE_POS,MATMUL_POS,INV_SCALE_POS,INV_SCALE_UVAL, + SCALE_ALL,MATMUL_ALL,INV_SCALE_ALL C C cd LOGICAL GEOM_SAVE INTEGER IANALYSE_REST,IWRITE_FIRST REAL ALPHA,ZERO REAL F_VALUE,STEP REAL FG0,FX0 LOGICAL LREFIN_SAVE,OCCUP_REF_FLAG_SAVE LOGICAL HCALC_SAVE,GRADC_SAVE LOGICAL CONV_FLAG,GAMMA_FLAG,STEP_FLAG LOGICAL LINMIN_FLAG CHARACTER LINE*128,MON_STYLE_SAVE*4 DATA ZERO/0.0/ C C---Initialisation step LINMIN_FLAG = .FALSE. GAMMA_FLAG = .FALSE. STEP_FLAG = .FALSE. GAMMA = 0.00 STEP = 0.05 IF(NTLS_CYCLE.GT.0.OR.TLS_FILE_FLAG) LINMIN_FLAG = .TRUE. HCALC = .TRUE. XCALC_FLAG = .TRUE. L_REMAIN = QQDEN-MAX(NVTMP,NVPOS) IERROR = 0 IF(L_REMAIN.LE.MAXATOM) THEN CALL ERRWRT(0,'Not enough memory. Change QQDEN') IERROR = 1 RETURN ENDIF MON_STYLE_SAVE = MON_STYLE OCCUP_REF_FLAG_SAVE = OCCUP_REF_FLAG DO NCYCL = 1,NCCREF IF(OCCUP_REF_FLAG_SAVE.AND.MOD(NCYCL,NCYCLE_OCC).EQ.0) THEN OCCUP_REF_FLAG = .TRUE. ELSE OCCUP_REF_FLAG = .FALSE. ENDIF IF((NCYCL.EQ.1.OR.NCYCL.EQ.NCCREF)) THEN MON_STYLE = MON_STYLE_SAVE ELSEIF(MON_STYLE.NE.'MANY') THEN MON_STYLE = 'FEW' ENDIF CALL ERRWRT(-1,' ') WRITE(LINE,'(5x,A,I6)')'CGMAT cycle number = ',NCYCL CALL ERRWRT(-1,' ') CALL ERRWRT(-1,LINE) IWRITE_FIRST = 1 IANALYSE_REST = 0 IF(REFID.NE.'UNRE') THEN IANALYSE_REST = 1 ENDIF ALPHA = 1.0 CALL ANALYSE_VDW(ALPHA,IWRITE_FIRST,IANALYSE_REST) NV1 = NVPOS + NVTMP CALL INIT_VEC(NV1,DV,ZERO) CALL RESIDVAL(F0) FG0 = FGEOM FX0 = FXRAY IF(MON_STYLE.NE.'NONE') THEN IF(GEOM_FLAG)CALL REPORT_OVER_GEOM IF(X_RAY_FLAG)CALL REPORT_XRAY_STATS ENDIF cd ALL_SCALES_FLAG = .FALSE. C C----Solution of equation with sparse matrix. IERROR = 0 C C---First solve equation for positional parameters NV1 = NVPOS TOLER = -1.0 NCONG = 500 cd IF(NDIS.EQ.0) THEN C C---No non-diagonls cd CALL DIAG_INVERSE_R(AM,V,DV) cd ELSE c IF(POS_REF_FLAG) THEN c WRITE(*,*)'In Pos ref' cd CALL CGSOLVE(NVPOS,QQV,NCONG,TOLER,AM(1),V(1),DV(1),VAL(1), cd + VAL(NVPOS+1),VAL(2*NVPOS+1),VAL(3*NVPOS+1), cd + VAL(4*NVPOS+1),L_REMAIN,SCALE_POS,MATMUL_POS, cd + SCALE_SHIFTS_ALL,INV_SCALE_POS,IERROR) NV1 = NVPOS+NVTMP cd CALL CGSOLVE(NV1,QQV,NCONG,TOLER,AM(1),V(1),DV(1),VAL(1), cd + VAL(NV1+1),VAL(2*NV1+1),VAL(3*NV1+1), cd + VAL(4*NV1+1),L_REMAIN,SCALE_ALL,MATMUL_ALL, cd + SCALE_SHIFTS_ALL,INV_SCALE_ALL,IERROR) C C---Condition close atoms C CALL CLOSE_ATOM_CONDITION(AM,VAL,QQDEN) C cd gamma = 0.01 cd gamma_flag = .true. cd DO IA=1,N_ATOM cd if(ia.le.20) write(*,*)atom_ref_flag(ia) cd IF(ATOM_REF_FLAG(IA).GT.0) THEN cd I = ATOM_REF_FLAG(IA)/10 cd WRITE(21,*)(AM(6*I-6+J),J=1,6) cd write(21,*) OCCUP(I), (V(3*I-3+k),K=1,3) cd ENDIF cd enddo C C---Solve equation for a while using direct invertion of "block" diagonals C cd call solve_unrestrained(dv) C CALL CGSOL_RM(NV1,QQV,NCONG,TOLER,AM(1),V(1),DV(1),VAL(1), + VAL(NV1+1),VAL(2*NV1+1),VAL(3*NV1+1), + VAL(4*NV1+1),L_REMAIN,SCALE_ALL,MATMUL_ALL, + SCALE_SHIFTS_ALL,INV_SCALE_ALL,IERROR,GAMMA_FLAG, + GAMMA,STEP_FLAG,STEP) cd stop C cd IF(GAMMA.GT.1.0) THEN cd GAMMA = 0.0 cd STEP = 0.05 cd NV1 = NVPOS+NVTMP cd NV = NVPOS cd CALL CGSOL_RM(NV,QQV,NCONG,TOLER,AM(1),V(1),DV(1),VAL(1), cd + VAL(NV+1),VAL(2*NV+1),VAL(3*NV+1), cd + VAL(4*NV+1),L_REMAIN,SCALE_POS,MATMUL_POS, cd + SCALE_SHIFTS_ALL,INV_SCALE_POS,IERROR,GAMMA_FLAG, cd + GAMMA,STEP_FLAG,STEP) cd GAMMA = 0.0 cd STEP = 0.05 cd NV = NVTMP cd WRITE(*,*)'U Values' cd CALL CGSOL_RM(NV,QQV,NCONG,TOLER,AM(NMPOS+9*NDIS+1),V(NVPOS+1), cd & DV(NVPOS+1), cd & VAL(1), cd + VAL(NV+1),VAL(2*NV+1),VAL(3*NV+1), cd + VAL(4*NV+1),L_REMAIN,SCALE_UVAL,MATMUL_UVAL, cd + SCALE_SHIFTS_ALL,INV_SCALE_UVAL,IERROR,GAMMA_FLAG, cd + GAMMA,STEP_FLAG,STEP) cd STOP cd DO I=1,10 cd WRITE(*,*)DV(I),XYZ_CRD(1,I),XYZ_CRD(2,I),XYZ_CRD(3,I) cd ENDDO cd WRITE(*,*) cd DO I=1,10 cd WRITE(*,*)DV(NVPOS+I), cd & U_ANISO(1,I),U_ANISO(2,I),U_ANISO(3,I), cd & U_ANISO(4,I),U_ANISO(5,I),U_ANISO(6,I) cd ENDDO C NV1 = NVPOS+NVTMP C cd STOP cd GAMMA = 0.0 cd STEP = 0.05 cd ENDIF IF(IERROR.NE.0) THEN CALL ERRWRT(0,'In CGSPARCE after CGSOLVE after positional'// + ' refinement') IERROR = 1 RETURN ENDIF c ENDIF C C---Now solve equations if we have B/U value refinements. c IF(ITEMP.NE.0.AND.REFID.NE.'IDEA') THEN c NVMAX = QQV-NVTMP c TOLER = -1.0 c NCONG = 500 c L_REMAIN = QQDEN - 5*NVTMP c NMTMP_POS = NMPOS+9*NDIS+1 c NVTMP_POS = NVPOS+1 c NMMAXU = QQM-NMTMP_POS+1 c NVMAXU = QQV-NVTMP_POS+1 c CALL CGSOLVE(NVTMP,NVMAX,NCONG,TOLER,AM(NMTMP_POS), c + V(NVTMP_POS),DV(NVTMP_POS),VAL(1),VAL(NVTMP+1), c + VAL(2*NVTMP+1),VAL(3*NVTMP+1),VAL(4*NVTMP+1), c + L_REMAIN,SCALE_UVAL,MATMUL_UVAL,SCALE_SHIFTS_ALL, c + INV_SCALE_UVAL,IERROR) IF(IERROR.NE.0) THEN CALL ERRWRT(0,'In CGSPARCE. CGSOLVE failed.'// + ' After thermal parameter refinement') IERROR = 2 RETURN ENDIF c ENDIF cd ENDIF C C---Occupancy refinement here cd CALL LINMINX(NO,1,NAME,F0,F1,AL,ALB) C C---Line minimisation should be added here. C C---Apply shifts using solution of normal equations F1 = F0 AL = 1.0 ALB = 1.0 ALC = 1.0 C C---Refine occupancy if necessary IF(OCCUP_REF_FLAG) THEN CALL OCCUPANCY_REFINE ENDIF C---------------------------------------- CALL SAVE_ALL_PARAMS('S') CALL ADD_SCALE_CONTR CALL APPLY(AL,ALB) CALL MAKE_U_POSITIVE IF(.NOT.LINMIN_FLAG)GOTO 100 UOVER_ATOM =0.0 C AL = 1.0 ALB = 1.0 HCALC_SAVE = HCALC HCALC = .FALSE. MON_STYLE_SAVE = MON_STYLE MON_STYLE = 'NONE' GRADC_SAVE = GRADC GRADC = .FALSE. DO I=1,3 CALL ANALYSE_VDW(ALPHA,IWRITE_FIRST,IANALYSE_REST) CALL RESIDVAL(F_VALUE) cd WRITE(*,*)F_VALUE,F0 IF(F_VALUE.GT.F0*1.0005) THEN AL = AL/2.0 ALB = ALB/2.0 CALL SAVE_ALL_PARAMS('R') CALL ADD_SCALE_CONTR CALL APPLY(AL,ALB) CALL MAKE_U_POSITIVE UOVER_ATOM = 0.0 GOTO 10 ENDIF GOTO 20 10 CONTINUE ENDDO 20 CONTINUE CALL SAVE_OVERALL_PARAMS('R') MON_STYLE = MON_STYLE_SAVE HCALC = HCALC_SAVE GRADC = GRADC_SAVE cd ALL_SCALES_FLAG = .FALSE. CALL CONVERGED_HKON(AL,ALB,CONV_FLAG) IF(CONV_FLAG) GOTO 30 100 CONTINUE ENDDO 30 CONTINUE write(*,*) IF(NLPRGO.NE.0) THEN C---Calculate geometry also cd GEOM_SAVE = GEOM_FLAG cd GEOM_FLAG = .FALSE. LREFIN_SAVE = LREFIN LREFIN = .FALSE. XCALC_FLAG = .TRUE. HCALC = .FALSE. GRADC = .FALSE. IWRITE_FIRST = 1 IANALYSE_REST = 0 IF(REFID.NE.'UNRE') THEN IANALYSE_REST = 1 ENDIF CALL ANALYSE_VDW(ALPHA,IWRITE_FIRST,IANALYSE_REST) CALL RESIDVAL(F1) IF(MON_STYLE.NE.'NONE') THEN IF(GEOM_FLAG )CALL REPORT_OVER_GEOM IF(X_RAY_FLAG)CALL REPORT_XRAY_STATS ENDIF cd GEOM_FLAG = GEOM_SAVE LREFIN = LREFIN_SAVE ENDIF OCCUP_REF_FLAG = OCCUP_REF_FLAG_SAVE C RETURN END C subroutine solve_unrestrained(dv) implicit none include 'atom_com.fh' include 'vitals.fh' include 'pls_incl.fh' include 'refi_flags.fh' real dv(*) REAL GX(3*MAXATOM),GQ(MAXATOM),GU1(6*MAXATOM) COMMON /REF_SPG/ GX,GU1,GQ C REAL H_XXD,H_UUD(3),H_UQD,H_QQD integer iv,ia,ierror character at_full1*12 C open(20) iv = 0 do ia=1,n_atom if(atom_ref_flag(ia).gt.0) then CALL FAST_HESSIAN_DIAGONAL(IA,H_XXD,H_UUD,H_UQD,H_QQD) dv(iv+1) = gx(iv+1)/h_xxd dv(iv+2) = gx(iv+2)/h_xxd dv(iv+3) = gx(iv+3)/h_xxd write(20,*)gx(iv+1),gx(iv+2),gx(iv+3) if(abs(dv(iv+1)).gt.1.0.or.abs(dv(iv+2)).gt.1.0.or. & abs(dv(iv+3)).gt.1.0) then call full_atom_name(ia,at_full1,ierror) write(*,*)at_full1,gx(iv+1),gx(iv+2),gx(iv+3) cd return endif iv = iv + 3 endif enddo cd write(*,*)iv cd write(*,*)dv(iv+1),dv(iv+2),dv(iv+3) C cd stop return end C C--- SUBROUTINE SAVE_ALL_PARAMS(SAVE_RESTORE) IMPLICIT NONE INCLUDE 'atom_com.fh' INCLUDE 'agreem.fh' INCLUDE 'weights.fh' INCLUDE 'refi_flags.fh' CHARACTER*1 SAVE_RESTORE C C---This subroutine saves all parameters for future use. For example C---in linear minimisation it may be necessary C INCLUDE 'save_all_params.fh' INTEGER IX,IA,IP,I C IF(SAVE_RESTORE.EQ.'S') THEN DO IA = 1,N_ATOM DO IX=1,3 XYZ_SAVE(IX,IA) = XYZ_CRD(IX,IA) ENDDO DO IX=1,6 U_ANISO_SAVE(IX,IA) = U_ANISO(IX,IA) ENDDO ENDDO DO IX=1,NBIN_ML SIGMA_ML_SAVE(IX) = SIGMA_ML(IX) DO IP=1,NPART+1 SCALE_ML_SAVE(IX,IP) = SCALE_ML(IX,IP) ENDDO ENDDO C D_ML_SCALE_OVER_SAVE = D_ML_SCALE_OVER D_ML_B_OVER_SAVE = D_ML_B_OVER D_ML_SCALE_BBULK_SAVE = D_ML_SCALE_BBULK D_ML_B_BBULK_SAVE = D_ML_B_BBULK DO I=1,NMAXP1 D_ML_SCALE_PART_SAVE(I) = D_ML_SCALE_PART(I) D_ML_B_PART_SAVE(I) = D_ML_B_PART(I) ENDDO DO IX=1,6 B_LS_ANISO_SAVE(IX) = B_LS_ANISO_OVER(IX) ENDDO SCALE_LS_OVER_SAVE = SCALE_LS_OVER B_LS_OVER_SAVE = B_LS_OVER SCALE_LS_BULK_SAVE = SCALE_LS_BULK B_LS_BULK_SAVE = B_LS_BULK DO IP=1,NPART SCALE_LS_PART_SAVE(IP) = SCALE_LS_PART(IP) B_LS_PART_SAVE(IP) = B_LS_PART(IP) ENDDO UOVER_SAVE = UOVER_ATOM C C--Save scale and ML parameters ELSE IF(SAVE_RESTORE.EQ.'R') THEN DO IA=1,N_ATOM DO IX = 1,3 XYZ_CRD(IX,IA) = XYZ_SAVE(IX,IA) ENDDO DO IX=1,6 U_ANISO(IX,IA) = U_ANISO_SAVE(IX,IA) ENDDO ENDDO DO IX=1,NBIN_ML SIGMA_ML(IX) = SIGMA_ML_SAVE(IX) DO IP=1,NPART+1 SCALE_ML(IX,IP) = SCALE_ML_SAVE(IX,IP) ENDDO ENDDO D_ML_SCALE_OVER = D_ML_SCALE_OVER_SAVE D_ML_B_OVER = D_ML_B_OVER_SAVE D_ML_SCALE_BBULK = D_ML_SCALE_BBULK_SAVE D_ML_B_BBULK = D_ML_B_BBULK_SAVE DO I=1,NMAXP1 D_ML_SCALE_PART(I) = D_ML_SCALE_PART_SAVE(I) D_ML_B_PART(I) = D_ML_B_PART_SAVE(I) ENDDO DO IX=1,6 B_LS_ANISO_OVER(IX) = B_LS_ANISO_SAVE(IX) ENDDO SCALE_LS_OVER = SCALE_LS_OVER_SAVE B_LS_OVER = B_LS_OVER_SAVE SCALE_LS_BULK = SCALE_LS_BULK_SAVE B_LS_BULK = B_LS_BULK_SAVE DO IP=1,NPART SCALE_LS_PART(IP) = SCALE_LS_PART_SAVE(IP) B_LS_PART(IP) = B_LS_PART_SAVE(IP) ENDDO UOVER_ATOM = UOVER_SAVE C C--Restore remaining parameters ELSE CALL ERRWRT(1,'Wrong argument for SAVE_ALL_PARAMETERS') ENDIF RETURN END C C--- SUBROUTINE SAVE_OVERALL_PARAMS(SAVE_RESTORE) IMPLICIT NONE INCLUDE 'atom_com.fh' INCLUDE 'agreem.fh' INCLUDE 'weights.fh' INCLUDE 'refi_flags.fh' CHARACTER*1 SAVE_RESTORE C C---This subroutine saves overall parameters for future use. For example C---in linear minimisation it may be necessary C INCLUDE 'save_all_params.fh' INTEGER IX,IA,IP,I C IF(SAVE_RESTORE.EQ.'S') THEN DO IX=1,NBIN_ML SIGMA_ML_SAVE(IX) = SIGMA_ML(IX) DO IP=1,NPART+1 SCALE_ML_SAVE(IX,IP) = SCALE_ML(IX,IP) ENDDO ENDDO D_ML_SCALE_OVER_SAVE = D_ML_SCALE_OVER D_ML_B_OVER_SAVE = D_ML_B_OVER D_ML_SCALE_BBULK_SAVE = D_ML_SCALE_BBULK D_ML_B_BBULK_SAVE = D_ML_B_BBULK DO I=1,NMAXP1 D_ML_SCALE_PART_SAVE(I) = D_ML_SCALE_PART(I) D_ML_B_PART_SAVE(I) = D_ML_B_PART(I) ENDDO DO IX=1,6 B_LS_ANISO_SAVE(IX) = B_LS_ANISO_OVER(IX) ENDDO SCALE_LS_OVER_SAVE = SCALE_LS_OVER B_LS_OVER_SAVE = B_LS_OVER SCALE_LS_BULK_SAVE = SCALE_LS_BULK B_LS_BULK_SAVE = B_LS_BULK DO IP=1,NPART SCALE_LS_PART_SAVE(IP) = SCALE_LS_PART(IP) B_LS_PART_SAVE(IP) = B_LS_PART(IP) ENDDO UOVER_SAVE = UOVER_ATOM C C--Save scale and ML parameters ELSE IF(SAVE_RESTORE.EQ.'R') THEN DO IX=1,NBIN_ML SIGMA_ML(IX) = SIGMA_ML_SAVE(IX) DO IP=1,NPART+1 SCALE_ML(IX,IP) = SCALE_ML_SAVE(IX,IP) ENDDO ENDDO D_ML_SCALE_OVER = D_ML_SCALE_OVER_SAVE D_ML_B_OVER = D_ML_B_OVER_SAVE D_ML_SCALE_BBULK = D_ML_SCALE_BBULK_SAVE D_ML_B_BBULK = D_ML_B_BBULK_SAVE DO I=1,NMAXP1 D_ML_SCALE_PART(I) = D_ML_SCALE_PART_SAVE(I) D_ML_B_PART(I) = D_ML_B_PART_SAVE(I) ENDDO DO IX=1,6 B_LS_ANISO_OVER(IX) = B_LS_ANISO_SAVE(IX) ENDDO SCALE_LS_OVER = SCALE_LS_OVER_SAVE B_LS_OVER = B_LS_OVER_SAVE SCALE_LS_BULK = SCALE_LS_BULK_SAVE B_LS_BULK = B_LS_BULK_SAVE DO IP=1,NPART SCALE_LS_PART(IP) = SCALE_LS_PART_SAVE(IP) B_LS_PART(IP) = B_LS_PART_SAVE(IP) ENDDO UOVER_ATOM = UOVER_SAVE C C--Restore remaining parameters ELSE CALL ERRWRT(1,'Wrong argument for SAVE_ALL_PARAMETERS') ENDIF RETURN END C SUBROUTINE CDIR c c---Minimization of functional by conjugate gradient method C C---Tronrud, Acta Cryst 1992 A48 pp 912-916 INCLUDE 'atom_com.fh' INCLUDE 'pls_incl.fh' INCLUDE 'cif_incl.fh' INCLUDE 'refi_flags.fh' INCLUDE 'agreem.fh' INCLUDE 'weights.fh' INCLUDE 'monitor.fh' INCLUDE 'const.fh' INCLUDE 'vitals.fh' c INTEGER NSYM_DIST(4,QQD) COMMON /LSQMTX / AM( QQM),V( QQV),DV( QQV) COMMON /DISTNS/ NE( QQD),NZZ( QQD),NWTB( QQD),NSYM_DIST C CHARACTER LINE*80 C C---Mtz PARAMETER (MCOLS = 200) INTEGER LOOKUP(MCOLS) COMMON /MTZRD/NLPRGI,NLPRGO,LOOKUP C LOGICAL GEOM_SAVE LOGICAL LREFIN_SAVE CALL ERRWRT(-1,'This option is dead') CALL ERRWRT(-1,'Use CGMAT instead') CALL ERRWRT(1,'Too little time to do everything') RETURN END C SUBROUTINE CGRAD c c---Minimization of functional by conjugate gradient method C C---See W.H.Press, B.P.Flannery, S.A.Teukolsky, W.T.Vetterling C--- Numerical recipes. The Art of Scientific Compputing. C--- pp 301-306 INCLUDE 'atom_com.fh' INCLUDE 'pls_incl.fh' INCLUDE 'cif_incl.fh' INCLUDE 'refi_flags.fh' INCLUDE 'agreem.fh' INCLUDE 'weights.fh' INCLUDE 'monitor.fh' INCLUDE 'const.fh' INCLUDE 'vitals.fh' INCLUDE 'ncs_rest.fh' INTEGER NSYM_DIST(4,QQD) COMMON /LSQMTX/ AM( QQM),V( QQV),DV( QQV) COMMON /DISTNS/ NE( QQD),NZZ( QQD),NWTB( QQD),NSYM_DIST C CHARACTER LINE*80 LOGICAL GEOM_SAVE LOGICAL LREFIN_SAVE C CALL ERRWRT(-1,'This option is dead') CALL ERRWRT(-1,'Use CGMAT instead') CALL ERRWRT(1,'Too little time to do everything') RETURN END C FUNCTION ESTIMF(AL1,CNST,IREF) C C---Estimate function value at AL*DV INCLUDE 'monitor.fh' INCLUDE 'refi_flags.fh' INCLUDE 'agreem.fh' C CHARACTER IREF*(*) LOGICAL SAVE_LOGIC C BMAX_SAVE = Bresetmax Bresetmax = 1.0E10 BMIN_SAVE = Bresetmin Bresetmin = -1.0E10 SAVE_LOGIC = GRADC GRADC = .FALSE. C IF(IREF.EQ.'XYZ') THEN AL = AL1 ALB = CNST ELSEIF(IREF.EQ.'BRF') THEN AL = CNST ALB = AL1 ELSE AL = AL1 ALB = CNST ENDIF C CALL APPLY(AL,ALB) CALL RESIDVAL(F) ESTIMF = F C C---Restore old parameters CALL APPLY(-AL,-ALB) FXRAY = FXSAVE GRADC = SAVE_LOGIC BresetMIN = BMIN_SAVE BresetMAX = BMAX_SAVE RETURN END C SUBROUTINE RESIDVAL(F) C C----Calculate value of residual (X_ray + geom) INCLUDE 'atom_com.fh' INCLUDE 'pls_incl.fh' INCLUDE 'weights.fh' INCLUDE 'refi_flags.fh' INCLUDE 'agreem.fh' INCLUDE 'cif_incl.fh' INCLUDE 'vitals.fh' INCLUDE 'monitor.fh' INCLUDE 'tls.fh' INCLUDE 'rharvest.fh' INCLUDE 'ncs_rest.fh' INTEGER NSYM_DIST(4,QQD) COMMON /LSQMTX/ AM( QQM),V( QQV),DV( QQV) COMMON /DISTNS/ NE( QQD),NZZ( QQD),NWTB( QQD),NSYM_DIST REAL GX(3*MAXATOM),GQ(MAXATOM),GU1(6*MAXATOM) COMMON /REF_SPG/ GX,GU1,GQ COMMON /R_SCRATCH/ VAL(QQDEN) INTEGER IT,NV1 REAL ALPHA,TLS_TRACE LOGICAL LINITIAL DATA ZERO/0.0/,IZERO/0/ C C---Restore vital parameters cd CALL RESTORE_VIT c c---Initialise vectors NNM = MIN(QQM,27*N_ATOM + 45*NDIS) CALL INIT_VEC(NNM,AM,ZERO) NV1 = NVPOS + NVTMP CALL INIT_VEC(NV1,V,ZERO) CALL INIT_VECN(NDIS,NE,IZERO) CALL INIT_VECN(NDIS,NZZ,IZERO) CALL INIT_VECN(NDIS,NWTB,IZERO) DO I=1,NDIS DO J=1,4 NSYM_DIST(J,I) = 0 ENDDO ENDDO cd CALL CIFINIT NRESTR = 0 C c---Make list of restraint for second derivatives. C---Call geom and refall/sfref for calculation C---geometric and x_ray terms respectively CALL MAKE_DIST_LIST F = 0.0 FGEOM = 0.0 IF(GEOM_FLAG) THEN HNRESTR = 0 CALL GEOM F = FGEOM ENDIF C C---No X-ray. We can return now. IF(.NOT.X_RAY_FLAG) RETURN DO IA = 1,N_ATOM DO I = 1,6 UANISO_OLD(I,IA) = U_ANISO(I,IA) ENDDO ENDDO IF(NTLS_CYCLE.GT.0.OR.TLS_FILE_FLAG) THEN LINITIAL = .FALSE. CALL TLS2ANISOU(LINITIAL) CALL MAKE_U_POSITIVE UOVER_ATOM = 0.0 ENDIF CALL REFALL if(nccref.le.0) return cd CALL ADD_TLS2DERIVS cd IF(NTLS_CYCLE.GT.0.OR.TLS_FILE_FLAG) cd & CALL ADD_TLS2DERIVS(GX,GU1,HX,HUU) C c---Add contibution of TLS here. WX = 1.0 WG = 1.0 C C---Ncs constraints here. C cd SUM_SAME = 0.0 cd CALL SAME_ATOMS(SUM_SAME) cd F = F + SUM_SAME CALL XGEOMADD(FXRAY,FGEOM,F,WX,WG,AM,V,GX,GU1,GQ, + VAL,QQDEN,IERROR) DO IA=1,N_ATOM DO I=1,6 U_ANISO(I,IA) = UANISO_OLD(I,IA) ENDDO ENDDO C C---Treat close atoms. I.e. if atoms are too close put restraints so that C---they will try to become same C IF(NUMBER_NCSC.GT.0) CALL NCS_CONST(QQM,QQV,AM,V,VAL,QQDEN) c ENDIF C C---Convert derivatives of U values to derivatives of S, where C---U = S^T S + alpha I. S is triangular matrix. This reparameterisation C---ensures that U is positive definite. When applying shifts this fact C---should be remembered. This reparametrisation might change behaviour of C---function cd CALL CONVERT_DERIVS_TO_S C RETURN END C SUBROUTINE CLOSE_ATOM_CONDITION(AM,U_HIST,NWORKSPACE) IMPLICIT NONE C C--- Check and try to stabilise cases when atoms are too close c--- each other C INCLUDE 'atom_com.fh' INCLUDE 'pls_incl.fh' INCLUDE 'celsym.fh' INCLUDE 'const.fh' INCLUDE 'vitals.fh' INCLUDE 'refi_flags.fh' INCLUDE 'agreem.fh' INCLUDE 'monitor.fh' INCLUDE 'weights.fh' INCLUDE 'tls.fh' C REAL AM(*) INTEGER NWORKSPACE INTEGER U_HIST(*) C INTEGER N_OBJECT(QQD),N_TARGET(QQD),NW_UVAL(QQD),NSYM_DIST(4,QQD) COMMON /DISTNS/ N_OBJECT,N_TARGET,NW_UVAL,NSYM_DIST C C--- local variables INTEGER I,J,K INTEGER ID,IA1,IA2,IA INTEGER IPOS,IA11 INTEGER COUNT_UP REAL CLOSE_DIST,DIST1 INTEGER LWORK,INFO REAL AMAT(3,3),AMAT1(3,3),E_V(3),WORKSPACE(40) REAL AMAX_AM,AM_ADD REAL AISO_ADD REAL CALC_DISTANCE EXTERNAL CALC_DISTANCE C LWORK = 40 INFO = 0 CLOSE_DIST = 0.6 AISO_ADD = 800.0 C C---Find positions of U values in the diagonal of hessian COUNT_UP = NMPOS + 9*NDIS IA1 = 0 DO IA=1,N_ATOM IF(ATOM_REF_FLAG(IA).GT.0) THEN IA1 = IA1 + 1 U_HIST(IA1) = COUNT_UP IF(U_ANISO(2,IA).LE.0.0) THEN COUNT_UP = COUNT_UP + 1 ELSE COUNT_UP = COUNT_UP + 21 ENDIF ENDIF ENDDO C C---loop over distances DO ID=1,NDIS IA1 = N_OBJECT(ID) IA2 = N_TARGET(ID) DIST1 = CALC_DISTANCE(IA1,IA2,NSYM_DIST(1,ID)) IF(DIST1.LE.CLOSE_DIST) THEN C C--- Distance between atoms is too close. Sort them out C IF(IA1.EQ.IA2) THEN C C-- Same atoms. Only one of them IA11 = ATOM_REF_FLAG(IA1)/10 IF(IA11.GT.0) THEN IPOS = 6*IA11-5 cd AMAX_AM = AMAX1(AM(IPOS), cd & AMAX1(AM(IPOS+1),AM(IPOS+2))) C C---Take care if all of them are 0 or small cd WRITE(*,*) AMAX_AM cd AM_ADD = AMAX_AM*(CLOSE_DIST-DIST1)/CLOSE_DIST AMAT(1,1) = AM(IPOS) AMAT(2,2) = AM(IPOS+1) AMAT(3,3) = AM(IPOS+2) AMAT(1,2) = AM(IPOS+3) AMAT(1,3) = AM(IPOS+4) AMAT(2,3) = AM(IPOS+5) AMAT(2,1) = AMAT(1,2) AMAT(3,1) = AMAT(1,3) AMAT(3,2) = AMAT(2,3) cd WRITE(*,*) AMAT CALL SSYEV('V','U',3,AMAT,3,E_V,WORKSPACE,LWORK,INFO) E_V(1) = AMAX1(E_V(1),E_V(3)*(1.0-DIST1/CLOSE_DIST)*0.8) E_V(2) = AMAX1(E_V(2),E_V(3)*(1.0-DIST1/CLOSE_DIST)*0.8) DO I=1,3 DO J=1,3 AMAT1(I,J) = 0.0 DO K=1,3 AMAT1(I,J) = AMAT1(I,J) + & AMAT(K,I)*AMAT(K,J)*E_V(K) ENDDO ENDDO ENDDO cd AM(IPOS ) = AM(IPOS ) + AM_ADD cd AM(IPOS+1) = AM(IPOS+1) + AM_ADD cd AM(IPOS+2) = AM(IPOS+2) + AM_ADD AM(IPOS ) = AMAT1(1,1) AM(IPOS+1) = AMAT1(2,2) AM(IPOS+2) = AMAT1(3,3) AM(IPOS+3) = AMAT1(1,2) AM(IPOS+4) = AMAT1(1,3) AM(IPOS+5) = AMAT1(2,3) C C---U values IF(U_ANISO(2,IA1).LE.0.0) THEN C--- c--- Do nothing. Isotropic case cd WRITE(*,*)AM(U_HIST(IA1)) C---Isotropic case. cd IF(AM(U_HIST(IA1)).GT.100.0) AISO_ADD = AM(U_HIST(IA1)) cd AM(U_HIST(IA1)) = AM(U_HIST(IA1)) + cd & (CLOSE_DIST-DIST1)/CLOSE_DIST*AISO_ADD ELSE C C---Anisotropic case. Eigenvaluse should be analysed. IPOS = U_HIST(IA1) AMAX_AM = AISO_ADD AMAX_AM = -1.0E32 DO I=1,6 AMAX_AM = AMAX1(AMAX_AM,AM(IPOS+I)) ENDDO AM_ADD = AMAX_AM*(CLOSE_DIST-DIST1)/CLOSE_DIST DO I=1,6 AM(IPOS+I) = AM(IPOS+I) + AM_ADD ENDDO ENDIF ENDIF ELSE GOTO 100 C C--- Different atoms. Sort out diagonal terms only. IA11 = ATOM_REF_FLAG(IA1)/10 IF(IA11.GT.0) THEN IPOS = 6*IA11-5 AMAX_AM = AMAX1(AM(IPOS), & AMAX1(AM(IPOS+1),AM(IPOS+2))) cd WRITE(*,*) AMAX_AM AM_ADD = AMAX_AM*(CLOSE_DIST-DIST1)/CLOSE_DIST AM(IPOS ) = AM(IPOS ) + AM_ADD AM(IPOS+1) = AM(IPOS+1) + AM_ADD AM(IPOS+2) = AM(IPOS+2) + AM_ADD C C---U values IF(U_ANISO(2,IA1).LE.0.0) THEN C---Isotropic case. do nothing. IF(AM(U_HIST(IA1)).GT.100.0) AISO_ADD = AM(U_HIST(IA1)) AM(U_HIST(IA1)) = AM(U_HIST(IA1)) + & (CLOSE_DIST-DIST1)/CLOSE_DIST*AISO_ADD ELSE C C---Anisotropic case. IPOS = U_HIST(IA1) cd AMAX_AM = AISO_ADD AMAX_AM = -1.0E32 DO I=1,6 AMAX_AM = AMAX1(AMAX_AM,AM(IPOS+I)) ENDDO AM_ADD = AMAX_AM*(CLOSE_DIST-DIST1)/CLOSE_DIST DO I=1,6 AM(IPOS+I) = AM(IPOS+I) + AM_ADD ENDDO ENDIF ENDIF C C---Second atom IA11 = ATOM_REF_FLAG(IA2)/10 IF(IA11.GT.0) THEN IPOS = 6*IA11-5 AMAX_AM = AMAX1(AM(IPOS), & AMAX1(AM(IPOS+1),AM(IPOS+2))) cd WRITE(*,*) AMAX_AM AM_ADD = AMAX_AM*(CLOSE_DIST-DIST1)/CLOSE_DIST AM(IPOS ) = AM(IPOS ) + AM_ADD AM(IPOS+1) = AM(IPOS+1) + AM_ADD AM(IPOS+2) = AM(IPOS+2) + AM_ADD C C---U values IF(U_ANISO(2,IA1).LE.0.0) THEN C---Isotropic case. This position might be empty. So process it. IF(AM(U_HIST(IA2)).GT.100.0) AISO_ADD = AM(U_HIST(IA2)) AM(U_HIST(IA2)) = AM(U_HIST(IA2)) + & (CLOSE_DIST-DIST1)/CLOSE_DIST*AISO_ADD ELSE C C---Anisotropic case. IPOS = U_HIST(IA2) cd AMAX_AM = AISO_ADD AMAX_AM = -1.0E32 DO I=1,6 AMAX_AM = AMAX1(AMAX_AM,AM(IPOS+I)) ENDDO AM_ADD = AMAX_AM*(CLOSE_DIST-DIST1)/CLOSE_DIST DO I=1,6 AM(IPOS+I) = AM(IPOS+I) + AM_ADD ENDDO ENDIF ENDIF ENDIF C C---U values C ENDIF 100 CONTINUE ENDDO RETURN END C C REAL FUNCTION CALC_DISTANCE(IA1,IA2,NSYM_DIST) IMPLICIT NONE C INCLUDE 'atom_com.fh' INCLUDE 'celsym.fh' C INTEGER I,IX INTEGER IS,ITX(3) INTEGER IA1,IA2,NSYM_DIST(4) REAL XYZ_TMP(3),XYZ_TMP1(3) REAL XYZ1(3),XYZ2(3) REAL DIST1 LOGICAL ERROR C IS = NSYM_DIST(1) IF(IS.LE.0) THEN IS = 1 ITX(1) = 0 ITX(2) = 0 ITX(3) = 0 ELSE ITX(1) = NSYM_DIST(2) ITX(2) = NSYM_DIST(3) ITX(3) = NSYM_DIST(4) ENDIF IF(IS.GT.1.OR.ITX(1).NE.0.OR.ITX(2).NE.0.OR.ITX(3).NE.0) THEN CALL MAT2VEC(3,3,CS_ORT_TO_FRAC,XYZ_CRD(1,IA2),XYZ_TMP,ERROR) DO I=1,3 XYZ_TMP1(I) = RealSymmMatrx(I,1,IS)*XYZ_TMP(1) + & RealSymmMatrx(I,2,IS)*XYZ_TMP(2) + & RealSymmMatrx(I,3,IS)*XYZ_TMP(3) + & RealSymmMatrx(I,4,IS) + FLOAT(ITX(I)) ENDDO CALL MAT2VEC(3,3,CS_FRAC_TO_ORT,XYZ_TMP1,XYZ2,ERROR) ELSE DO IX=1,3 XYZ2(IX) = XYZ_CRD(IX,IA2) ENDDO ENDIF C DO IX=1,3 XYZ1(IX) = XYZ_CRD(IX,IA1) ENDDO DIST1 = SQRT(ABS((XYZ1(1)-XYZ2(1))**2+(XYZ1(2)-XYZ2(2))**2+ & (XYZ1(3)-XYZ2(3))**2)) CALC_DISTANCE = DIST1 END SUBROUTINE ADD_TLS2DERIVS(GX,GU1) C C---Adds contribution of TLS to individual atomic derivatives. At the C---moment only the fact that atoms could be refined isotropicly but c---due to TLS they could appear anisotropicly so derivatives for anisotropic C---atoms will be calculated and they have to be converted to isotropic c---case. In more general case correlation between positional and thermal C---parameters should be taken into account also (More precisely as U values C---are dependent on positional parameters U = U(x), derivatives of C---positional parameters will be affected by U) IMPLICIT NONE REAL GX(*),GU1(*) INCLUDE 'atom_com.fh' INCLUDE 'tls.fh' C INTEGER IA,I,IPOSV1,IPOSV2,IA1,IA11 REAL DU_ISO,D2U_ISO,HUU_MAX C IPOSV1 = 0 IPOSV2 = 0 HUU_MAX = -1.0E32 DO IA=1,N_ATOM IF(ATOM_REF_FLAG(IA).GT.0) THEN IA1 = ATOM_REF_FLAG(IA)/10 IA11 = ATOM_REF_FLAG(IA)-IA1*10 IF(IA11.GE.3) THEN IF(UANISO_OLD(2,IA).NE.0.0) THEN C c----Atom is anisotropic copy gradients DO I=1,6 IPOSV1 = IPOSV1 + 1 IPOSV2 = IPOSV2 + 1 GU1(IPOSV1) = GU1(IPOSV2) ENDDO ELSEIF(UANISO_OLD(2,IA).LE.0.0.AND.IGROUP(IA).GT.0) THEN C c---Atom is isotropic DU_ISO = & (GU1(IPOSV2+1)+GU1(IPOSV2+2)+GU1(IPOSV2+3)) IPOSV2 = IPOSV2 + 6 IPOSV1 = IPOSV1 + 1 GU1(IPOSV1) = DU_ISO ELSEIF(UANISO_OLD(2,IA).LE.0.0) THEN IPOSV1 = IPOSV1 + 1 IPOSV2 = IPOSV2 + 1 GU1(IPOSV1) = GU1(IPOSV2) ENDIF ELSEIF(UANISO_OLD(2,IA).LE.0.0.AND.IGROUP(IA).EQ.0) THEN IPOSV1 = IPOSV1 + 1 IPOSV2 = IPOSV2 + 1 GU1(IPOSV1) = 0.0 ELSEIF(UANISO_OLD(2,IA).EQ.0.0.AND.IGROUP(IA).NE.0) THEN IPOSV1 = IPOSV1 + 1 IPOSV2 = IPOSV2 + 6 GU1(IPOSV1) = 0.0 ELSEIF(UANISO_OLD(2,IA).NE.0.0) THEN DO I=1,6 IPOSV1 = IPOSV1 + 1 IPOSV2 = IPOSV2 + 1 GU1(IPOSV1) = 0.0 ENDDO ENDIF ENDIF ENDDO C RETURN END C SUBROUTINE ADD_TLS2DERIVS1(GX,GU1,HX,HUU) C C---Adds contribution of TLS to individual atomic derivatives. At the C---moment only the fact that atoms could be refined isotropicly but c---due to TLS they could appear anisotropicly so derivatives for anisotropic C---atoms will be calculated and they have to be converted to isotropic c---case. In more general case correlation between positional and thermal C---parameters should be taken into account also (More precisely as U values C---are dependent on positional parameters U = U(x), derivatives of C---positional parameters will be affected by U) IMPLICIT NONE REAL GX(*),GU1(*),HX(*),HUU(*) INCLUDE 'atom_com.fh' INCLUDE 'tls.fh' C INTEGER IA,I,IPOSM1,IPOSM2,IPOSV1,IPOSV2,IA1,IA11 REAL DU_ISO,D2U_ISO,HUU_MAX C cd DO IA=1,N_ATOM cd DO I=1,6 cd U_ANISO(I,IA) = UANISO_OLD(I,IA) cd ENDDO cd ENDDO IPOSM1 = 0 IPOSM2 = 0 IPOSV1 = 0 IPOSV2 = 0 HUU_MAX = -1.0E32 DO IA=1,N_ATOM IF(ATOM_REF_FLAG(IA).GT.0) THEN IA1 = ATOM_REF_FLAG(IA)/10 IA11 = ATOM_REF_FLAG(IA)-IA1*10 IF(IA11.GE.3) THEN IF(UANISO_OLD(2,IA).GT.0.0) THEN C c----Atom is anisotropic copy gradients and second derivatives DO I=1,6 IPOSV1 = IPOSV1 + 1 IPOSV2 = IPOSV2 + 1 GU1(IPOSV1) = GU1(IPOSV2) ENDDO DO I=1,21 IPOSM1 = IPOSM1 + 1 IPOSM2 = IPOSM2 + 1 HUU(IPOSM1) = HUU(IPOSM2) ENDDO ELSEIF(UANISO_OLD(2,IA).LE.0.0.AND.IGROUP(IA).GT.0) THEN C c---Atom is isotropic DU_ISO = & (GU1(IPOSV2+1)+GU1(IPOSV2+2)+GU1(IPOSV2+3)) IPOSV2 = IPOSV2 + 6 IPOSV1 = IPOSV1 + 1 GU1(IPOSV1) = DU_ISO D2U_ISO = (HUU(IPOSM2+1)+HUU(IPOSM2+2)+HUU(IPOSM2+3) & +2.0*(HUU(IPOSM2+7)+HUU(IPOSM2+8)+HUU(IPOSM2+12))) IPOSM1 = IPOSM1 + 1 HUU(IPOSM1) = D2U_ISO IPOSM2 = IPOSM2 + 21 ELSEIF(UANISO_OLD(2,IA).LE.0.0) THEN IPOSV1 = IPOSV1 + 1 IPOSV2 = IPOSV2 + 1 IPOSM1 = IPOSM1 + 1 IPOSM2 = IPOSM2 + 1 GU1(IPOSV1) = GU1(IPOSV2) HUU(IPOSM1) = HUU(IPOSM2) ENDIF ELSEIF(UANISO_OLD(2,IA).EQ.0.0.AND.IGROUP(IA).EQ.0) THEN IPOSV1 = IPOSV1 + 1 IPOSV2 = IPOSV2 + 1 GU1(IPOSV1) = 0.0 IPOSM1 = IPOSM1 + 1 IPOSM2 = IPOSM2 + 1 HUU(IPOSM1)= 0.0 ELSEIF(UANISO_OLD(2,IA).EQ.0.0.AND.IGROUP(IA).NE.0) THEN IPOSV1 = IPOSV1 + 1 IPOSV2 = IPOSV2 + 6 GU1(IPOSV1) = 0.0 IPOSM1 = IPOSM1 + 1 IPOSM2 = IPOSM2 + 21 HUU(IPOSM1)= 0.0 ELSEIF(UANISO_OLD(2,IA).NE.0.0) THEN DO I=1,6 IPOSV1 = IPOSV1 + 1 IPOSV2 = IPOSV2 + 1 GU1(IPOSV1) = 0.0 ENDDO DO I=1,21 IPOSM1 = IPOSM1 + 1 IPOSM2 = IPOSM2 + 1 HUU(IPOSM1) = 0.0 ENDDO ENDIF ENDIF ENDDO C RETURN END C SUBROUTINE FINAL_CHECK_XRAY_DERIVS(HX,HU) IMPLICIT NONE INCLUDE 'atom_com.fh' INCLUDE 'refi_flags.fh' INCLUDE 'vitals.fh' REAL HX(*),HU(*) C C---Tries to make derivatives positive. It might happen with current version INTEGER IA1,IA11,IA,IU,IPAR,IPARU REAL HX_MAX,HU_MAX,HU_MAX1,HX_MAX1 C HX_MAX = -1.0E32 HU_MAX = -1.0E32 IPAR = 0 IPARU = 0 DO IA=1,N_ATOM IF(ATOM_REF_FLAG(IA).GT.0) THEN IA1 = ATOM_REF_FLAG(IA)/10 IA11 = ATOM_REF_FLAG(IA)- IA1*10 IF(IA11.GT.2) THEN HX_MAX = AMAX1(HX_MAX,AMAX1(HX(IPAR+1), & AMAX1(HX(IPAR+2),HX(IPAR+3)))) IPAR = IPAR + 6 IF(ITEMP.GT.0) THEN IF(U_ANISO(2,IA).EQ.0.0) THEN HU_MAX = AMAX1(HU_MAX,HU(IPARU+1)) IPARU = IPARU + 1 ELSE HU_MAX = AMAX1(HU_MAX,AMAX1(HU(IPARU+1), & AMAX1(HU(IPARU+2), & AMAX1(HU(IPARU+3),AMAX1(HU(IPARU+4), & AMAX1(HU(IPARU+5),HU(IPARU+6))))))) IPARU = IPARU + 21 ENDIF ENDIF ELSE IPAR = IPAR + 6 IF(ITEMP.GT.0) THEN IF(U_ANISO(2,IA).EQ.0.0) THEN IPARU = IPARU + 1 ELSE IPARU = IPARU + 21 ENDIF ENDIF ENDIF ENDIF ENDDO C HX_MAX1 = HX_MAX HU_MAX1 = HU_MAX IPAR = 0 IPARU = 0 DO IA=1,N_ATOM IF(ATOM_REF_FLAG(IA).GT.0) THEN IA1 = ATOM_REF_FLAG(IA)/10 IA11 = ATOM_REF_FLAG(IA)- IA1*10 IF(IA11.GT.2) THEN IF(HX(IPAR+1).LE.0.0.OR.HX(IPAR+2).LE.0.0.OR. & HX(IPAR+3).LE.0.0) THEN HX(IPAR+1) = HX_MAX1 HX(IPAR+2) = HX_MAX1 HX(IPAR+3) = HX_MAX1 HX(IPAR+4) = 0.0 HX(IPAR+5) = 0.0 HX(IPAR+6) = 0.0 ENDIF IPAR = IPAR + 6 IF(ITEMP.GT.0) THEN IF(U_ANISO(2,IA).EQ.0.0) THEN IF(HU(IPARU+1).LE.0.0) HU(IPARU+1) = HU_MAX1 IPARU = IPARU + 1 ELSE IF(HU(IPARU+1).LE.0.0.OR.HU(IPARU+2).LE.0.0.OR. & HU(IPARU+4).LE.0.0.OR.HU(IPARU+5).LE.0.0.OR. & HU(IPARU+6).LE.0.0) THEN DO IU=1,6 HU(IPARU+IU) = HU_MAX1 ENDDO DO IU=7,21 HU(IPARU+IU) = 0.0 ENDDO ENDIF IPARU = IPARU + 21 ENDIF ENDIF ELSE IPAR = IPAR + 6 IF(ITEMP.GT.0) THEN IF(U_ANISO(2,IA).EQ.0.0) THEN IPARU = IPARU + 1 ELSE IPARU = IPARU + 21 ENDIF ENDIF ENDIF ENDIF ENDDO RETURN END C--- SUBROUTINE ERRWRT(I,LINE) C C---SUBROUTINE writes message line to output file an C---I=0 - warning message does not exit c---otherwise Error message exits CHARACTER LINE*(*) COMMON /FILES / ISYSR,ISYSW,IFOFC,ISCRF,IATMR,IDISK,JDISK,NEWFIL . ,INXYZ,IOUTX,MTZIN C L = MAX(1,LENSTR(LINE)) IF(I.EQ.-1) THEN WRITE(ISYSW,'(A)')LINE(1:L) ELSEIF(I.EQ.0) THEN WRITE(ISYSW,'(A,A)')'===> Warning: ',LINE(1:L) RETURN ELSEIF(I.EQ.1) THEN WRITE(ISYSW,'(A,A)')'===> Error: ',LINE(1:L) CALL REFMAC_CLEAN_UP_FILES CALL CCPERR(1,LINE) ENDIF END C SUBROUTINE APPLY(STEP,STEPB) IMPLICIT NONE c c---Subroutine applies shifts to parameters INCLUDE 'atom_com.fh' INCLUDE 'pls_incl.fh' INCLUDE 'celsym.fh' INCLUDE 'celsym_aniso.fh' INCLUDE 'const.fh' INCLUDE 'refi_flags.fh' INCLUDE 'agreem.fh' INCLUDE 'monitor.fh' INCLUDE 'vitals.fh' INCLUDE 'occupancy_params.fh' REAL STEP,STEPB REAL AM,V,DV COMMON /LSQMTX/ AM( QQM),V( QQV),DV( QQV) INTEGER IA,IANISO,LV,IVPOS real c1,c2 real shift_allow_high,shift_allow_low real rms_shift,ave_shift,max_shift REAL DVX,DVY,DVZ,DBDF REAL U_DELTA(6) REAL U_MIN,UEIGEN REAL UANISO0(3,3),UANISO1(6),EVALUE(3),WORKSPACE(20) INTEGER LWORK,INFO,I,II,IERROR integer inegatives,irefine real shift_correct C LWORK = 20 C C---First analyse if U shifts are too large. If so then reduce shifts C inegatives = 0 irefine = 0 shift_correct = 1.0 IF(ITEMP.GT.0) THEN c2 = BDAMP*STEPB LV = NVPOS + 1 DO IA=1,N_ATOM IF(ATOM_REF_FLAG(IA).GT.0) THEN irefine = irefine + 1 IF(U_ANISO(2,IA).LE.0.0) THEN C C--- Isotropic U-values DBDF = DV(LV)*c2 IF(.NOT.(EXCLUDE_OCC_B.AND.OCCUP_REF_FLAG)) THEN IF(U_ANISO(1,IA) + DBDF.LE.BresetMin) then inegatives = inegatives + 1 cd WRITE(*,*)U_ANISO(1,IA)+DBDF ENDIF LV = LV + 1 ENDIF ELSE C C---Anisotropic U-values DO IANISO=1,6 U_DELTA(IANISO) = DV(LV) LV = LV + 1 ENDDO IF(.NOT.(EXCLUDE_OCC_B.AND.OCCUP_REF_FLAG)) THEN DO IANISO=1,6 DBDF = c2*U_DELTA(IANISO) UANISO1(IANISO) = U_ANISO(IANISO,IA) + DBDF ENDDO C C---Make U values positive UANISO0(1,1) = UANISO1(1) UANISO0(2,2) = UANISO1(2) UANISO0(3,3) = UANISO1(3) UANISO0(1,2) = UANISO1(4) UANISO0(1,3) = UANISO1(5) UANISO0(2,3) = UANISO1(6) UANISO0(2,1) = UANISO0(1,2) UANISO0(3,1) = UANISO0(1,3) UANISO0(3,2) = UANISO0(2,3) CALL SSYEV('V','U',3,UANISO0,3,EVALUE,WORKSPACE, & LWORK,INFO) IF(EVALUE(1).LT.BresetMin.OR. & EVALUE(2).LT.BresetMin.OR. & EVALUE(3).LT.BresetMin) THEN inegatives = inegatives + 1 cd WRITE(*,*)EVALUE(1),EVALUE(2),EVALUE(3) ENDIF ENDIF ENDIF ENDIF ENDDO if(inegatives.gt.1) shift_correct = 0.5 cd WRITE(*,*)inegatives,inegatives/irefine ENDIF cd stop C C-- Positional C-- C C--- First find maximum, average and rms shifts. They can be used to modify C-- shift lengths C c1 = STEP*PDAMP*shift_correct IF(NVPOS.GT.0) THEN call find_avemaxrms(NVPOS,DV,ave_shift,max_shift,rms_shift) cd write(*,*)'average, maximum and rms ', cd & ave_shift*c1,max_shift*c1,rms_shift*c1 C C---Now use this information for shift readjustment IVPOS = 0 shift_allow_high = (ave_shift + 10.0*rms_shift)*c1 shift_allow_low = (ave_shift - 10.0*rms_shift)*c1 shift_allow_high = amax1(0.0,amin1(shift_allow_high,0.35)) shift_allow_low = amin1(0.0,amax1(shift_allow_low,-0.35)) shift_allow_high = 0.5 shift_allow_low = -0.5 DO IA=1,N_ATOM IF(ATOM_REF_FLAG(IA).GT.0) THEN DVX = DV(IVPOS+1)*c1 if(DVX.gt.shift_allow_high) DVX = shift_allow_high if(DVX.lt.shift_allow_low ) DVX = shift_allow_low DVY = DV(IVPOS+2)*c1 if(DVY.gt.shift_allow_high) DVY = shift_allow_high if(DVY.lt.shift_allow_low ) DVY = shift_allow_low DVZ = DV(IVPOS+3)*c1 if(DVZ.gt.shift_allow_high) DVZ = shift_allow_high if(DVZ.lt.shift_allow_low ) DVZ = shift_allow_low IVPOS = IVPOS + 3 XYZ_CRD(1,IA) = XYZ_CRD(1,IA) + DVX XYZ_CRD(2,IA) = XYZ_CRD(2,IA) + DVY XYZ_CRD(3,IA) = XYZ_CRD(3,IA) + DVZ ENDIF ENDDO ENDIF C C--- Thermal C--- Similarly for B values. Shifts should not bee too large C--- c2 = BDAMP*STEPB*shift_correct c1 = c2*PISQ8 IF(ITEMP.GT.0) THEN C C--- Find average, maximum and rms shifts. How to define tese terms for C--- anisotropic ADPs C call find_avemaxrms(NVTMP,DV(NVPOS+1),ave_shift,max_shift, & rms_shift) cd write(*,*)'average, maximum and rms B shifts ', cd * ave_shift*c1,max_shift*c1,rms_shift*c1 C C--- Now use this inofrmation for shift readjustment C LV = NVPOS + 1 shift_allow_high = ( ave_shift + 10.0*rms_shift)*c2 shift_allow_low = ( ave_shift - 10.0*rms_shift)*c2 shift_allow_high = 20.0/PISQ8 shift_allow_low = -20.0/PISQ8 DO IA=1,N_ATOM IF(ATOM_REF_FLAG(IA).GT.0) THEN IF(U_ANISO(2,IA).LE.0.0) THEN C C--- Isotropic U-values DBDF = DV(LV)*c2 IF(.NOT.(EXCLUDE_OCC_B.AND.OCCUP_REF_FLAG)) THEN if(DBDF.gt.shift_allow_high) DBDF = shift_allow_high if(DBDF.lt.shift_allow_low ) DBDF = shift_allow_low U_ANISO(1,IA) = U_ANISO(1,IA) + DBDF U_ANISO(1,IA) = AMAX1(BresetMin, * AMIN1(U_ANISO(1,IA),BResetMax)) ENDIF LV = LV + 1 ELSE C C---Anisotropic U-values DO IANISO=1,6 U_DELTA(IANISO) = DV(LV) LV = LV + 1 ENDDO IF(.NOT.(EXCLUDE_OCC_B.AND.OCCUP_REF_FLAG)) THEN DO IANISO=1,6 DBDF = c2*U_DELTA(IANISO) if(DBDF.gt.shift_allow_high) DBDF = shift_allow_high if(DBDF.lt.shift_allow_low ) DBDF = shift_allow_low U_ANISO(IANISO,IA) = U_ANISO(IANISO,IA) + DBDF ENDDO C C---Make U values positive UANISO0(1,1) = U_ANISO(1,IA) UANISO0(2,2) = U_ANISO(2,IA) UANISO0(3,3) = U_ANISO(3,IA) UANISO0(1,2) = U_ANISO(4,IA) UANISO0(1,3) = U_ANISO(5,IA) UANISO0(2,3) = U_ANISO(6,IA) UANISO0(2,1) = UANISO0(1,2) UANISO0(3,1) = UANISO0(1,3) UANISO0(3,2) = UANISO0(2,3) CALL SSYEV('V','U',3,UANISO0,3,EVALUE,WORKSPACE, & LWORK,INFO) EVALUE(1) = AMAX1(BResetMin,AMIN1(BResetMax,EVALUE(1))) EVALUE(2) = AMAX1(BResetMin,AMIN1(BResetMax,EVALUE(2))) EVALUE(3) = AMAX1(BResetMin,AMIN1(BResetMax,EVALUE(3))) CALL EIGEN2U(EVALUE,UANISO0,UANISO1) DO IANISO=1,6 U_ANISO(IANISO,IA) = UANISO1(IANISO) ENDDO ENDIF ENDIF ENDIF ENDDO ENDIF C cd CALL CHECK_U_VALUES cd U_MIN = 1.E32 cd DO IA=1,N_ATOM cd IF(ATOM_REF_FLAG(IA).GT.0) THEN cd IF(U_ANISO(2,IA).LE.0.0) THEN cd U_MIN = AMIN1(U_MIN,U_ANISO(1,IA)) cd ELSE cd CALL FIND_MIN_EIGEN(U_ANISO(1,IA),UEIGEN) cd U_MIN = AMIN1(U_MIN,UEIGEN) cd ENDIF cd ENDIF cd ENDDO cd WRITE(*,*)U_MIN,BResetMin cd IF(U_MIN.LE.0.0) THEN cd CALL ERRWRT(1,'It should not happend. negtive U') cd ENDIF C C---Add occupancy shifts. This routine does not know C---style of occupancy refineement IF(OCCUP_REF_FLAG) THEN LV = NVPOS + NVTMP + 1 DO IA=1,N_ATOM IF(ATOM_REF_FLAG(IA).GT.0) THEN OCCUP(IA) = AMAX1(0.0,AMIN1(1.0,OCCUP(IA)+QDAMP*DV(LV))) LV = LV + 1 ENDIF ENDDO ENDIF C C--- cd CALL CHECK_OCCUPANCIES RETURN END C SUBROUTINE CONVERGED_HKON(AL,ALB,CONV_FLAG) IMPLICIT NONE C INCLUDE 'refi_flags.fh' INCLUDE 'atom_com.fh' INCLUDE 'pls_incl.fh' INCLUDE 'vitals.fh' INCLUDE 'const.fh' C REAL AL,ALB LOGICAL CONV_FLAG C REAL AM,V,DV COMMON /LSQMTX/ AM( QQM),V( QQV),DV( QQV) INTEGER IA,IANISO,LV,IVPOS C REAL DPMAX,DBMAX C C--Positional c--This parameter should be replaced IF(NVPOS.GT.0) THEN IVPOS = 0 DPMAX = -1.0E32 DO IA=1,N_ATOM IF(ATOM_REF_FLAG(IA).GT.0) THEN DPMAX = AMAX1(DPMAX,ABS(DV(IVPOS+1)*AL)) DPMAX = AMAX1(DPMAX,ABS(DV(IVPOS+2)*AL)) DPMAX = AMAX1(DPMAX,ABS(DV(IVPOS+3)*AL)) IVPOS = IVPOS + 3 ENDIF ENDDO ENDIF cd WRITE(*,*)'MAximum positional shift ',DPMAX C C---Thermal IF(ITEMP.GT.0) THEN LV = NVPOS + 1 DBMAX = -1.0E32 DO IA=1,N_ATOM IF(ATOM_REF_FLAG(IA).GT.0) THEN IF(U_ANISO(2,IA).EQ.0.0) THEN C C---Isotropic U-values DBMAX = AMAX1(DBMAX,ABS(DV(LV)*ALB)) LV = LV + 1 ELSE C C---Anisotropic U-values DO IANISO=1,6 DBMAX = AMAX1(DBMAX,ABS(ALB*DV(LV))) LV = LV + 1 ENDDO ENDIF ENDIF ENDDO ENDIF cd WRITE(*,*)'Maximum B value shift ',DBMAX*PISQ8 C C--- CONV_FLAG = .FALSE. cd WRITE(*,*)'MAx shifts ',DPMAX,DBMAX IF(DPMAX.LE.EPS_CONV_X.AND.DBMAX.LE.EPS_CONV_B) THEN CONV_FLAG = .TRUE. ENDIF cd CALL CHECK_OCCUPANCIES RETURN END C SUBROUTINE SAVE_VITAL C C----Save vital parameters INCLUDE 'vitals.fh' C NAS = NA NDISS = NDIS NPLNS = NPLN NCHRS = NCHR NVDWS = NVDW NHBDS = NHBD NOCCS = NOCC NTORS = NTOR ITEMPS = ITEMP RETURN END C SUBROUTINE RESTORE_VIT C C---Restore vital parameters INCLUDE 'vitals.fh' C NA = NAS NDIS = NDISS NPLN = NPLNS NCHR = NCHRS NVDW = NVDWS NHBD = NHBDS NOCC = NOCCS NTOR = NTORS ITEMP = ITEMPS RETURN END C SUBROUTINE GEOM C C---Subroutine for calculation geometric part of C---residulas and its derivatives INCLUDE 'atom_com.fh' INCLUDE 'pls_incl.fh' INCLUDE 'celsym.fh' INCLUDE 'refi_flags.fh' INCLUDE 'agreem.fh' INCLUDE 'cif_incl.fh' INCLUDE 'const.fh' INCLUDE 'weights.fh' INCLUDE 'monitor.fh' INCLUDE 'vitals.fh' C C-----replace by big_mass.fh INTEGER NSYM_DIST(4,QQD) COMMON /LSQMTX/ AM( QQM),V( QQV),DV( QQV) COMMON /DISTNS/ N_OBJECT(QQD),N_TARGET(QQD),NW_UVAL(QQD), & NSYM_DIST C C---Local variables COMMON /R_SCRATCH/ WORKSPACE(QQDEN) C---Initialisation of first and second derivative matrices should be done C---before this routine NMTRX = QQM MXDIS = QQD C IF(PDEL.NE.0.0) THEN IMPOS = 0 DO IA=1,N_ATOM IF(ATOM_REF_FLAG(IA).GT.0) THEN AM(IMPOS+1) = 1.0/PDEL**2 AM(IMPOS+2) = 1.0/PDEL**2 AM(IMPOS+3) = 1.0/PDEL**2 IMPOS = IMPOS + 6 ENDIF ENDDO ENDIF c c---Contribution from geometry N_WORKSPACE = QQDEN-51*N_ATOM CALL CALC_GEOM_CONTR(WORKSPACE(1),WORKSPACE(N_ATOM+1), & WORKSPACE(51*N_ATOM+1),N_WORKSPACE) C RETURN END C SUBROUTINE HEADER(LABEL) CHARACTER LABEL*(*) CHARACTER*60 LINE C L = LEN(LABEL) LINE(1:4) = '****' LINE(57:60) = '****' L = 60 - 2*4 - LEN(LABEL) IF(L.LT.0) THEN CALL ERRWRT(0, ' TOO LONG CHARACTER STRING IN HEADER') RETURN ENDIF LB1 = L/2 LB2 = 4 + LB1 + LEN(LABEL) DO I=5,4+LB1 LINE(I:I)=' ' ENDDO DO I=LB2+1,56 LINE(I:I)=' ' ENDDO J=0 DO I=4+LB1+1,LB2 J=J+1 LINE(I:I)=LABEL(J:J) ENDDO CALL ERRWRT(-1,' ') CALL ERRWRT(-1,' '//LINE) CALL ERRWRT(-1,' ') RETURN END C SUBROUTINE REFALL C C--- IMPLICIT NONE INCLUDE 'atom_com.fh' INCLUDE 'solvent.fh' INCLUDE 'pls_incl.fh' C REAL VAL(QQDEN) COMMON /R_SCRATCH/VAL IF(SOLVENT_FLAG) CALL SOLVENT(VAL,QQDEN) CALL REFALL1(VAL,QQDEN) RETURN END C SUBROUTINE REFALL1(VAL,QQDEN) IMPLICIT NONE C----------------------------------------------------------------------- C THIS ROUTINE DESIGNED TO EMULATE SIMILAR ROUTINES FROM THE AGARWAL C FFT UNCONSTRAINED REFINEMENT PACKAGE. IT DRIVES THE SEQUENCE OF C FOURIER AND GRADIENT ROUTINES FOR THE CALCULATION OF LEAST-SQUARES C MATRIX ELEMENTS FOR PROLSQ. MODIFIED FOR ALL SPACE GROUPS C *********SPACE GROUP GENERAL********** C G.N.M. 16.09.91 C----------------------------------------------------------------------- INCLUDE 'celsym.fh' INCLUDE 'refi_flags.fh' INCLUDE 'agreem.fh' INCLUDE 'const.fh' INCLUDE 'vitals.fh' INCLUDE 'atom_com.fh' INCLUDE 'monitor.fh' INTEGER QQDEN REAL VAL(*) INTEGER HMAX,KMAX,LMAX COMMON /HKLLIM/ HMAX, KMAX, LMAX C REAL ASYMLIM1,ASYMLIM2,ASYMLIM3 INTEGER FC_ADDR,PHASE_ADDR,DEN_ADDR,IND_ADDR,FO_ADDR,SIGO_ADDR, + POOL_ADDR,FREE_ADDR,LIND_ADDR,NGX,NGY,NGZ,NMFOUR,NXY, & IREMADD,IESTIM,N_LAST,NWORKSPACE,NWORKSPACE1,NREF INTEGER SIZE,SIZEL CHARACTER LINE*128 C ----------------------------------------------------- C C---Find minimum accetable grid spacing for ffts NX = 0 NY = 0 NZ = 0 CALL GET_GRID_SPACING(FSHANN,NGX,NGY,NGZ,HMAX,KMAX,LMAX) CALL ASYLIM(CS_NSPGR,IPX,IPY,IPZ,NMFOUR) ASYMLIM1 = 1.0/REAL(IPX) ASYMLIM2 = 1.0/REAL(IPY) ASYMLIM3 = 1.0/REAL(IPZ) IF(MON_STYLE.EQ.'MANY') THEN WRITE(LINE,'(A,3F5.2)')' Limits of asymmetric unit :', + ASYMLIM1,ASYMLIM2,ASYMLIM3 CALL ERRWRT(-1,LINE) WRITE(LINE,'(A,3I5)') ' Grid spacing to be used : ', + NX,NY,NZ CALL ERRWRT(-1,LINE) WRITE(LINE,'(A,3I5)') ' Maximuum H,K,L : ', + HMAX,KMAX,LMAX CALL ERRWRT(-1,LINE) WRITE(LINE,'(A,3I5)') ' Minimum acceptable grid spacing: ', + NGX,NGY,NGZ CALL ERRWRT(-1,LINE) ENDIF c IF(MOD(NX,NMFOUR).NE.0.OR.MOD(NY,NMFOUR).NE.0 + .OR.MOD(NZ,NMFOUR).NE.0) THEN CALL ERRWRT(1,' Change grid spacings ') ENDIF C C----Now allocate memory for S.F. Den and so on N1 = NX N2 = NY N3 = NZ/IPZ IF(IPZ.GT.1) N3 = N3 + 1 SIZE = N1*N2*N3+1 NXY = (NX+1)*(NY+1) C C ALLOCATE MEMORY FOR FO,SIGO and IND IND_ADDR = 1 FO_ADDR = IND_ADDR + NOBS + 1 SIGO_ADDR = FO_ADDR + NOBS + 1 FREE_ADDR = SIGO_ADDR + NOBS + 1 C C---If we are going to use PHI FCi or other informations they could be added C---here C ALLOCATE MEMORY FOR FC FC_ADDR = FREE_ADDR + NOBS + 1 PHASE_ADDR = FC_ADDR + NOBS*(NPART+1) + 1 LIND_ADDR = PHASE_ADDR + NOBS*(NPART+1) + 1 IF(MIR_FLAG.OR.PHASE_FLAG) LIND_ADDR = LIND_ADDR + 4*NOBS + 1 DEN_ADDR = PHASE_ADDR + NOBS + 1 POOL_ADDR = DEN_ADDR + SIZE + 1 C C---Check if size of intermediate arrays are within size of R_SCRATCH IREMADD = QQDEN-(POOL_ADDR + 2*NXY) IESTIM = POOL_ADDR + 2*NXY IF(IREMADD.LT.0) THEN CALL ERRWRT(-1,' ') CALL ERRWRT(-1,' ') WRITE(LINE,'(A,2I13)')' Size for R_SCRATCH is not enough '// + ' You need QQDEN at least ',IESTIM,QQDEN CALL ERRWRT(-1,LINE) CALL ERRWRT(1,'Increase QQDEN in pls_incl.fh and recompile ') ENDIF N_LAST = POOL_ADDR+2*NXY NWORKSPACE = QQDEN-POOL_ADDR NWORKSPACE1 = QQDEN - N_LAST CALL REF_ALL(NREF,VAL(IND_ADDR),VAL(FO_ADDR),VAL(SIGO_ADDR), * VAL(DEN_ADDR),VAL(FC_ADDR),VAL(PHASE_ADDR),VAL(FREE_ADDR), * VAL(POOL_ADDR),VAL(LIND_ADDR),NOBS, & VAL(N_LAST+1),NWORKSPACE,NWORKSPACE1) C RETURN END C C----------------------------------------------------------------------- SUBROUTINE REF_ALL(NREF,NIND,FO,SIGO,DEN,FC,PHASE,FREER, + WORKSPACE,LIND,SZ,WORKSPACE1,NWORKSPACE,NWORKSPACE1) IMPLICIT NONE C C THIS ROUTINE FROM THE AGARWAL FFT UNCONSTRAINED REFINEMENT C PACKAGE. MODIFIED FOR ALL SPACE GROUPS. C C----This subroutine uses semi FFT C C ********SPACE GROUP GENERAL*********** C G.N.M. 16.09.91 C----------------------------------------------------------------------- INCLUDE 'atom_com.fh' INCLUDE 'pls_incl.fh' INCLUDE 'celsym.fh' INCLUDE 'celsym_aniso.fh' INCLUDE 'weights.fh' INCLUDE 'refi_flags.fh' INCLUDE 'agreem.fh' INCLUDE 'vitals.fh' INCLUDE 'const.fh' INCLUDE 'hessian_impl.fh' INTEGER NREF INTEGER HMAX, KMAX, LMAX COMMON /HKLLIM/ HMAX, KMAX, LMAX INTEGER NWORKSPACE,NWORKSPACE1 INTEGER SZ,SZ1 REAL WORKSPACE(*),WORKSPACE1(*) REAL FC(*),PHASE(*),FO(*),SIGO(*),FREER(*),DEN(N1,N2,N3) REAL GX(3*MAXATOM),GQ(MAXATOM),GU1(6*MAXATOM) INTEGER NIND(*),LIND(*) COMMON /REF_SPG/ GX,GU1,GQ INTEGER IZERO,IPP,I,J,IQ,IA REAL ZERO integer ir,ihh(3) real umin_loc,umax_loc,u_add_loc,grid_max,umin_req real b_add_loc,rsq real lstlsq external lstlsq,unpack cd REAL HXX(6*MAXATOM),HUU(21*MAXATOM),HQQ(MAXATOM),HQU(6*MAXATOM) cd COMMON /REF_SEC/HXX,HUU,HQQ,HQU c INTEGER ierror,atom_1,atom_2,k,m,n,ID INTEGER N_PAIRS,COUNT,RATIO,JUMP cd REAL HW(9,50000),HBB(50000),HWMED cd COMMON /COMPARE/ HW,HBB cd COMMON /SUMERR/ SUMERRXX,SUMERRYY,SUMERRZZ,SUMERRXY,SUMERRXZ, cd & SUMERRYZ,SUMERRB cd REAL SUMERRXX(50000),SUMERRYY(50000),SUMERRZZ(50000), cd & SUMERRXY(50000),SUMERRXZ(50000),SUMERRYZ(50000), cd & SUMERRB(50000) REAL HPOSMED,DIST,BSUM,HBBD integer N_TARGET(QQD),N_OBJECT(QQD),NW_UVAL(QQD),NSYM_DIST(4,QQD) INTEGER VAL_VAR,START COMMON /DISTNS/ N_OBJECT,N_TARGET,NW_UVAL,NSYM_DIST character*12 atm_name_1,atm_name_2 REAL H_XXD,H_UUD(3),H_UQD,H_QQD REAL H_XX(3,3),H_UU(6,6),H_QQ,H_XB(3),H_XU(3,6),H_XQ(3),H_QX(3), & H_BQ,H_QB,H_UQ(6),H_QU(6),HWZ(9),HBBZ REAL AVXX, AVYY, AVZZ, AVXY, AVXZ, AVYZ, AVBB character at_full1*12 integer ia1,iv,i1,i2,i3 LOGICAL ERROR C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C----Initialise ZERO = 0.0 IZERO = 0 SZ1 = SZ*(NPART+1) CALL INIT_VEC(SZ1,FC,ZERO) CALL INIT_VEC(SZ1,PHASE,ZERO) CALL INIT_VEC(SZ,FO,ZERO) CALL INIT_VEC(SZ,SIGO,ZERO) CALL INIT_VECN(SZ,NIND,IZERO) CALL INIT_VECN(SZ,LIND,IZERO) C C GENERATE MODEL ELECTRON DENSITY MAP C C--Find minimum B values C CALL FIND_U_EXTREMES(umin_loc,umax_loc) grid_max = amax1(CS_CELL(1)/float(NX),amax1(CS_CELL(2)/float(NY), & CS_CELL(3)/float(NZ))) umin_req = grid_max**2/1.1 C C---If u values are too small then add some value to all atoms to avoid u_add_loc = 0.0 if(umin_req.gt.umin_loc) then u_add_loc = umin_req - umin_loc cd u_add_loc = 0.0 do ia=1,n_atom if(atom_ref_flag(ia).gt.0.0) then if(u_aniso(2,ia).le.0.0) then u_aniso(1,ia) = u_aniso(1,ia) + u_add_loc else u_aniso(1,ia) = u_aniso(1,ia) + u_add_loc u_aniso(2,ia) = u_aniso(2,ia) + u_add_loc u_aniso(3,ia) = u_aniso(3,ia) + u_add_loc endif endif enddo endif cd write(*,*)u_add_loc,umin_loc,umax_loc cd stop CALL DENSTY(DEN) C C Read refelctions (indices and observed s.f.) CALL RDIND(NREF,FO,SIGO,FREER,NIND) C C-----Reduce electron density to asymmetric unit (REDALLI) calculate C-----structure factors (RFFT) and read FO and SIGO and partial structure C-----factors (RDIND1) CALL REDALLI(N1,N2,N3,NX,NY,NZ,ROTR,TRR,NonCubSym,DEN) CALL RFFT(NREF,DEN,WORKSPACE,FC,PHASE,NIND) CALL RDIND1(NREF,FO,SIGO,FREER,FC,PHASE,NIND) C C--Now remove added b values from Fs b_add_loc = u_add_loc*pisq8 do ir=1,nref cd if(sigo(ir).gt.0.0) then CALL UNPACK(NIND(ir),ihh(1),ihh(2),ihh(3)) rsq = lstlsq(1,ihh(1),ihh(2),ihh(3)) fc(NPART*NOBS + ir) = fc(NPART*NOBS +ir)*exp(rsq*b_add_loc) cd endif enddo call write_reflections_r(nref,fo,sigo,fc,phase,nind) C C----Write structure factors to intermediate file C CALL WRTIND(NREF,FO,SIGO,FREER,PHASE,NIND) C C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C REPORT AGREEMENT FACTORS AND ACCUMULATE FOURIER COEFFICIENTS C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - CALL PROPI(NREF,NIND,FO,SIGO,FC,PHASE,FREER,WORKSPACE,NWORKSPACE) C C Last cycle - only calculate and output SFS IF( .NOT.LREFIN) THEN do ia=1,n_atom if(atom_ref_flag(ia).gt.0.0) then if(u_aniso(2,ia).le.0.0) then u_aniso(1,ia) = u_aniso(1,ia) - u_add_loc else u_aniso(1,ia) = u_aniso(1,ia) - u_add_loc u_aniso(2,ia) = u_aniso(2,ia) - u_add_loc u_aniso(3,ia) = u_aniso(3,ia) - u_add_loc endif endif enddo CALL MTZWRITE(NREF,LIND,NIND,FO,SIGO,FC,PHASE,FREER) RETURN ENDIF C C-----Gradient calculation using fft IF(GRADC) THEN IPP = 0 do ir=1,nref cd if(sigo(ir).gt.0.0) then CALL UNPACK(NIND(ir),ihh(1),ihh(2),ihh(3)) rsq = lstlsq(1,ihh(1),ihh(2),ihh(3)) fc(ir) = fc(ir)*exp(rsq*b_add_loc) cd sigo(ir) = sigo(ir)*exp(-rsq*b_add_loc) cd endif enddo C C---add b vlues for the coefficients for gradient also CALL FFT(NREF,DEN,FC,PHASE,NIND,IPP) CALL GRAD_ALL(DEN,GX,GU1,GQ) ENDIF C C---now remove b values from atoms C do ia=1,n_atom if(atom_ref_flag(ia).gt.0.0) then if(u_aniso(2,ia).le.0.0) then u_aniso(1,ia) = u_aniso(1,ia) - u_add_loc else u_aniso(1,ia) = u_aniso(1,ia) - u_add_loc u_aniso(2,ia) = u_aniso(2,ia) - u_add_loc u_aniso(3,ia) = u_aniso(3,ia) - u_add_loc endif endif enddo IF(HCALC) THEN cd CALL HDIAG_R(NREF,NIND,FO,HXX,HUU,HQQ,HQU) cd DO I=1,10 cd WRITE(*,*)(HXX(6*(I-1) + J),J=1,6),HUU(I) cd ENDDO CALL D2DA_RADIAL(NREF,NIND,FO,SIGO) cd CALL HDIAG_R(NREF,NIND,FO,HXX,HUU,HQQ,HQU) c cd print*, 'tabulation started' CALL FAST_HESSIAN_TABULATION(SMINB_ML(1),SMAXB_ML(NBIN_ML)) cd print*, 'tabulation ended' CALL CALCULATE_DPI_ML c ENDIF RETURN END C subroutine find_u_extremes(umin,umax) implicit none include 'atom_com.fh' include 'const.fh' real umin,umax integer ia,ig,il real u_iso,u_all_l C umax = -1.0E32 umin = 1.0E32 do ia=1,n_atom if(occup(ia).gt.0.0.and.atom_ref_flag(ia).gt.0) then il = id_sf(ia) if(u_aniso(2,ia) .le. 0.0) then u_iso = u_aniso(1,ia) else u_iso = (u_aniso(1,ia)+u_aniso(2,ia)+u_aniso(3,ia))/3.0 endif do ig=1,ngaus u_all_l = cs_b(ig,il)/PISQ8 + u_iso umin = amin1(umin,u_all_l) umax = amax1(umax,u_all_l) enddo endif enddo return end C SUBROUTINE CALCULATE_DPI_ML IMPLICIT NONE C C---This subroutine calculates DPIs (or approximate standard uncertainties) C---using likilohood function INCLUDE 'atom_com.fh' INCLUDE 'rharvest.fh' INCLUDE 'const.fh' C INTEGER IA,N_REFINED,IA1,IA11 REAL H_XXD,H_UUD(3),H_UQD,H_QQD HESU_ML = 0.0 HESU_ML_B = 0.0 N_REFINED = 0 DO IA=1,N_ATOM IA1 = ATOM_REF_FLAG(IA)/10 IA11 = ATOM_REF_FLAG(IA)-IA1*10 IF(IA11.GE.3) THEN N_REFINED = N_REFINED + 1 CALL FAST_HESSIAN_DIAGONAL(IA,H_XXD,H_UUD,H_UQD,H_QQD) HESU_ML = HESU_ML + H_XXD HESU_ML_B = HESU_ML_B + H_UUD(1) ENDIF ENDDO HESU_ML = SQRT(2.0*REAL(N_REFINED)/HESU_ML) HESU_ML_B = PISQ8*SQRT(2.0*REAL(N_REFINED)/HESU_ML_B) END C SUBROUTINE GET_GRID_SPACING(FSHANN,NGX,NGY,NGZ,HMAX,KMAX,LMAX) C c---Calculates grid spacing using CELL resolution,and shannon's rate C INCLUDE 'refi_flags.fh' INCLUDE 'celsym.fh' INCLUDE 'atom_com.fh' INTEGER HMAX,KMAX,LMAX SMAX = STLMAX*2 cd CALL LRCELL(1,CS_CELL) SHH = 1.0/CS_CELL(1)**2 SKK = 1.0/CS_CELL(2)**2 SLL = 1.0/CS_CELL(3)**2 HMAX = INT(CS_CELL(1)*SMAX) KMAX = INT(CS_CELL(2)*SMAX) LMAX = INT(CS_CELL(3)*SMAX) C C---SZ is size for Structure factors and so on RH0MAX = SMAX*CS_CELL(1) RK0MAX = SMAX*CS_CELL(2) RL0MAX = SMAX*CS_CELL(3) NGX = INT(RH0MAX+SQRT(.25/SHH+RH0MAX**2))+1 NGY = INT(RK0MAX+SQRT(.25/SKK+RK0MAX**2))+1 NGZ = INT(RL0MAX+SQRT(.25/SLL+RL0MAX**2))+1 C C----If number of grids were not given calculate them CALL ASYLIM(CS_NSPGR,IPX,IPY,IPZ,NMFOUR) IF(NX*NY*NZ.EQ.0) THEN NX = INT(NGX*FSHANN) NY = INT(NGY*FSHANN) NZ = INT(NGZ*FSHANN) C C---Now find best numbers CALL NXYZ235(NX,NMFOUR) CALL NXYZ235(NY,NMFOUR) CALL NXYZ235(NZ,NMFOUR) ELSE C C---Check if predefined number of grid points is less than minimum acceptable IBAD = 0 IF(NX.LT.NGX) IBAD = IBAD+1 IF(NY.LT.NGY) IBAD = IBAD+1 IF(NZ.LT.NGZ) IBAD = IBAD+1 IF(IBAD.GE.1) THEN NX = INT(NGX*FSHANN) NY = INT(NGY*FSHANN) NZ = INT(NGZ*FSHANN) C C---Now find best numbers CALL NXYZ235(NX,NMFOUR) CALL NXYZ235(NY,NMFOUR) CALL NXYZ235(NZ,NMFOUR) ENDIF ENDIF RETURN END C SUBROUTINE SET_HKON_FLAGS C IMPLICIT NONE INCLUDE 'atom_com.fh' INCLUDE 'refi_flags.fh' C C---Define atoms which should be included in refinement, in geometry C---calculation. In this version all atoms but hydrogens will be included c---in refinement. hydrogens will be included in geometry calculation c---and contribution of them will be taken into account in structure factor C---calculation INTEGER IA1,IA REAL U_ISOT C IA1 = 0 DO IA=1,N_ATOM IF(OCCUP(IA).LE.0.00001) THEN ATOM_REF_FLAG(IA) = 0 ELSEIF(CS_ELEMENT(ID_SF(IA)).EQ.'H ' + .OR.CS_ELEMENT(ID_SF(IA)).EQ.'H-1 ') THEN IA1 = IA1 + 1 IF(.NOT.HYDROGEN_REFINE_FLAG) THEN ATOM_REF_FLAG(IA) = IA1*10 + 2 ELSE IF(HYDROGEN_BVALUE.EQ.'ISOT') THEN ATOM_REF_FLAG(IA) = IA1*10 + 3 IF(U_ANISO(2,IA).NE.0.0) THEN U_ISOT = (U_ANISO(1,IA)+U_ANISO(2,IA)+U_ANISO(3,IA))/3.0 U_ANISO(1,IA) = U_ISOT U_ANISO(2,IA) = 0.0 U_ANISO(3,IA) = 0.0 U_ANISO(4,IA) = 0.0 U_ANISO(5,IA) = 0.0 U_ANISO(6,IA) = 0.0 ENDIF ELSE ATOM_REF_FLAG(IA) = IA1*10 + 4 ENDIF ENDIF ELSE IA1 = IA1 + 1 ATOM_REF_FLAG(IA) = IA1*10 + 4 ENDIF ENDDO C C---Change style of refinement for some of preselected atoms C C---No occupancy refinement DO IA=1,N_ATOM OCCUP_REF(IA) = 0 ENDDO C C---Some other checks and refinement settings should be added here RETURN END C SUBROUTINE AUTO_DEFINE_OCCUPANCY_GROUPS IMPLICIT NONE C INCLUDE 'atom_com.fh' INCLUDE 'occupancy_params.fh' INTEGER ICH,IRES,IATOM,IR_FIRST_C,IR_LAST_C, IA_FIRST_C,IA_LAST_C, & IDEFINE,NALTS,IGRP1,IGRP2,IALTS,IATOM1,NGROUP_OCC1, & NALTS_OCCUP_NEW(MAXALT_GROUPS), & NGROUP_ALT_NEW(MAXALTS,MAXALT_GROUPS),IGRP,NGROUP_OCC_OLD CHARACTER ALT_CODE_THIS C C---Initialise DO IATOM=1,N_ATOM IGRP_OCCUP(IATOM) = 0 ENDDO C C--First find refinable groups which have alt codes. NGROUP_OCC = 0 DO ICH=1,N_GROUP IR_FIRST_C = IRES_FIRST(ICH) IR_LAST_C = IRES_FIRST(ICH)+NRES_CHAIN(ICH)-1 DO IRES = IR_FIRST_C,IR_LAST_C IA_FIRST_C = IRATM_FIRST(IRES) IA_LAST_C = IA_FIRST_C + NATM_RES(IRES)-1 1 CONTINUE IDEFINE = 0 DO IATOM = IA_FIRST_C,IA_LAST_C IF(IGRP_OCCUP(IATOM).LE.0.AND. & ATOM_REF_FLAG(IATOM).GT.0) THEN IF(ID_ALT(IATOM).NE.'.') THEN NGROUP_OCC = NGROUP_OCC + 1 ALT_CODE_THIS = ID_ALT(IATOM) IGRP_OCCUP(IATOM)= NGROUP_OCC IDEFINE = 1 GOTO 20 ENDIF ENDIF ENDDO C 20 CONTINUE IF(IDEFINE.EQ.1) THEN DO IATOM=IA_FIRST_C,IA_LAST_C IF(IGRP_OCCUP(IATOM).LE.0.AND. & ATOM_REF_FLAG(IATOM).GT.0.AND. & ID_ALT(IATOM).EQ.ALT_CODE_THIS) THEN IGRP_OCCUP(IATOM) = NGROUP_OCC ENDIF ENDDO GOTO 1 ENDIF ENDDO ENDDO C C---Now define which group is alternative for which DO IGRP1=1,NGROUP_OCC NALTS_OCCUP(IGRP1) = 0 ENDDO DO ICH=1,N_GROUP IR_FIRST_C = IRES_FIRST(ICH) IR_LAST_C = IRES_FIRST(ICH)+NRES_CHAIN(ICH)-1 DO IRES = IR_FIRST_C,IR_LAST_C IA_FIRST_C = IRATM_FIRST(IRES) IA_LAST_C = IA_FIRST_C + NATM_RES(IRES)-1 100 CONTINUE IDEFINE = 0 DO IATOM = IA_FIRST_C,IA_LAST_C IF(IGRP_OCCUP(IATOM).GT.0) THEN IGRP1 = IGRP_OCCUP(IATOM) DO IATOM1=IA_FIRST_C,IA_LAST_C IF(IGRP_OCCUP(IATOM1).GT.0) THEN IF(ATM_NAME(IATOM).EQ.ATM_NAME(IATOM1).AND. & ID_ALT(IATOM).NE.ID_ALT(IATOM1)) THEN C C---Check if it is already in the list NALTS = NALTS_OCCUP(IGRP1) IGRP2 = IGRP_OCCUP(IATOM1) IF(NALTS.GT.0) THEN DO IALTS = 1,NALTS IF(IGRP2.EQ.NGROUP_ALT(IALTS,IGRP1)) GOTO 120 ENDDO ENDIF NALTS = NALTS + 1 NALTS_OCCUP(IGRP1) = NALTS NGROUP_ALT(NALTS,IGRP1) = IGRP2 120 CONTINUE ENDIF ENDIF ENDDO ENDIF ENDDO ENDDO ENDDO C C---Identify unique groups. DO IGRP = 1,NGROUP_OCC IF(NGROUP_OCC1.LE.0) THEN NGROUP_OCC1 = 1 NALTS_OCCUP_NEW(NGROUP_OCC1) = NALTS_OCCUP(IGRP) NGROUP_REF(NGROUP_OCC1) = IGRP DO IALTS=1,NALTS_OCCUP(IGRP) NGROUP_ALT_NEW(IALTS,NGROUP_OCC1) = NGROUP_ALT(IALTS,IGRP) ENDDO ELSE C C---Check if we already have this group DO IGRP1=1,NGROUP_OCC IF(IGRP.EQ.NGROUP_REF(IGRP1)) GOTO 300 DO IALTS=1,NALTS_OCCUP_NEW(IGRP1) IF(IGRP.EQ.NGROUP_ALT_NEW(IALTS,IGRP1)) GOTO 300 ENDDO ENDDO NGROUP_OCC1 = NGROUP_OCC1 + 1 NALTS_OCCUP_NEW(NGROUP_OCC1) = NALTS_OCCUP(IGRP) NGROUP_REF(NGROUP_OCC1) = IGRP IF(NALTS_OCCUP(IGRP).GT.0) THEN DO IALTS=1,NALTS_OCCUP(IGRP) NGROUP_ALT_NEW(IALTS,NGROUP_OCC1) = NGROUP_ALT(IALTS,IGRP) ENDDO ENDIF ENDIF 300 CONTINUE ENDDO NGROUP_OCC_OLD = NGROUP_OCC NGROUP_OCC = NGROUP_OCC1 DO IGRP = 1,NGROUP_OCC NALTS_OCCUP(IGRP) = NALTS_OCCUP_NEW(IGRP) DO IALTS=1,NALTS_OCCUP(IGRP) NGROUP_ALT(IALTS,IGRP) = NGROUP_ALT_NEW(IALTS,IGRP) ENDDO ENDDO C C---Now redefne groups using bond list. CALL BOND_2_OCCUPANCY(NGROUP_OCC_OLD) C C---Treat non bonded contacts and use them also C CALL NONBOND_2_OCCUPANCY cd STOP RETURN END C SUBROUTINE BOND_2_OCCUPANCY(NGROUP_OCC_OLD) IMPLICIT NONE INCLUDE 'atom_com.fh' INCLUDE 'occupancy_params.fh' INCLUDE 'restr_files.fh' INCLUDE 'pls_incl.fh' INTEGER N_TARGET,N_OBJECT,N_TEMPOR,NSYM_DIST(4,QQD) COMMON /DISTNS/ N_OBJECT( QQD),N_TARGET(QQD),N_TEMPOR(QQD), & NSYM_DIST INTEGER N_RELATE1(1000),N_RELATE2(1000),N_CLASSES(1000) C REAL RS_VIDL,RS_SDI INTEGER ISCRB,IFAIL,LL,IA1,IA2,ISYM1,ITX(3),N_BONDS,N_EQUIV, & IGRP1,IGRP2,IE,IB,NGROUP_OCC_OLD,IO,IO1,IO2,NALTS1, & NALTS2,IGP3,IGP4,N_EQUIV1,IE1,NGROUP_OCC_NEW,IGRP,NALTS INTEGER LENSTR EXTERNAL LENSTR ISCRB = 0 IFAIL = -1 IF(LENSTR(BOND_FILE).LE.0) RETURN CALL CCPDPN(ISCRB,BOND_FILE(1:LENSTR(BOND_FILE)),'UNKNOWN','U', & LL,IFAIL) C c--Read bond files and make list N_BONDS = 0 1 CONTINUE READ(ISCRB,END=99)IA1,IA2,RS_VIDL,RS_SDI,ISYM1,ITX(1),ITX(2), & ITX(3) IF(RS_SDI.LE.0) GOTO 1 IF(IA1.LE.0.OR.IA2.LE.0) GOTO 1 IF(ATOM_REF_FLAG(IA1).LE.0.OR.ATOM_REF_FLAG(IA2).LE.0) GOTO 1 cd IF(I_RESID(IA1).EQ.I_RESID(IA2)) GOTO 1 N_BONDS = N_BONDS + 1 N_TARGET(N_BONDS) = IA1 N_OBJECT(N_BONDS) = IA2 C C---Another bond? GOTO 1 99 CONTINUE CLOSE(ISCRB) C C--Now check bonds and make equivalence relations N_EQUIV = 0 DO IB=1,N_BONDS IA1 = N_TARGET(IB) IA2 = N_OBJECT(IB) IGRP1 = IGRP_OCCUP(IA1) IGRP2 = IGRP_OCCUP(IA2) IF(IGRP1.GT.0.AND.IGRP2.GT.0.AND.IGRP1.NE.IGRP2) THEN C c---Check if this relation is already in the list IF(N_EQUIV.GT.0) THEN DO IE = 1,N_EQUIV IF((N_RELATE1(IE).EQ.IGRP1.AND.N_RELATE2(IE).EQ.IGRP2).OR. & (N_RELATE1(IE).EQ.IGRP2.AND.N_RELATE2(IE).EQ.IGRP1)) & GOTO 199 ENDDO ENDIF C C--Not in the list WRITE(*,*)RES_NUM_PDB(I_RESID(IA1)),RES_NUM_PDB(I_RESID(IA2)), & IGRP1,IGRP2 N_EQUIV = N_EQUIV + 1 N_RELATE1(N_EQUIV) = IGRP1 N_RELATE2(N_EQUIV) = IGRP2 C 199 CONTINUE ENDIF ENDDO C C---Alternatives also. N_EQUIV1 = N_EQUIV DO IE=1,N_EQUIV IGRP1 = N_RELATE1(N_EQUIV) IGRP2 = N_RELATE2(N_EQUIV) DO IO = 1,NGROUP_OCC IF(NGROUP_REF(IO).EQ.IGRP1) THEN NALTS1 = NALTS_OCCUP(IO) IO1 = IO GOTO 200 ENDIF ENDDO 200 CONTINUE IF(NALTS1.LE.0) GOTO 300 DO IO=1,NGROUP_OCC IF(NGROUP_REF(IO).EQ.IGRP2) THEN NALTS2 = NALTS_OCCUP(IO) IO2 = IO GOTO 210 ENDIF ENDDO 210 CONTINUE IF(NALTS2.LE.0) GOTO 300 IF(NALTS1.EQ.1.AND.NALTS2.EQ.1) THEN IGP3 = NGROUP_ALT(1,IO1) IGP4 = NGROUP_ALT(1,IO2) DO IE1=1,N_EQUIV1 IF((N_RELATE1(IE1).EQ.IGP3.AND.N_RELATE2(IE1).EQ.IGP4).OR. & (N_RELATE1(IE1).EQ.IGP4.AND.N_RELATE2(IE1).EQ.IGP3)) & GOTO 300 ENDDO N_EQUIV1 = N_EQUIV1 + 1 N_RELATE1(N_EQUIV1) = IGP3 N_RELATE2(N_EQUIV1) = IGP4 C C---Different treatement if NALTS > 1 ENDIF 300 CONTINUE ENDDO N_EQUIV = N_EQUIV1 WRITE(*,*)'Number of equivalence relations',N_EQUIV,N_BONDS CALL EQUIVALENCE_CLASSES(NGROUP_OCC_OLD,N_EQUIV,N_RELATE1, & N_RELATE2, & N_TARGET) CALL TREE_2_CLASS_NUMBER(NGROUP_OCC_OLD,N_TARGET,N_CLASSES) DO IB=1,NGROUP_OCC_OLD WRITE(*,*)IB,N_CLASSES(IB) ENDDO C C--Now redefine groups DO IO=1,N_ATOM IF(IGRP_OCCUP(IO).GT.0) THEN IGRP_OCCUP(IO) = N_CLASSES(IGRP_OCCUP(IO)) ENDIF ENDDO NGROUP_REF(1) = N_CLASSES(NGROUP_REF(1)) NGROUP_OCC_NEW = 1 NALTS = NALTS_OCCUP(1) IF(NALTS.GT.0) THEN DO IE=1,NALTS NGROUP_ALT(IE,1) = N_CLASSES(NGROUP_ALT(IE,1)) ENDDO ENDIF DO IGRP=2,NGROUP_OCC IGRP1 = NGROUP_REF(IGRP) DO IGP3 = 1,NGROUP_OCC_NEW IF(NGROUP_REF(IGP3).EQ.N_CLASSES(IGRP1)) GOTO 400 ENDDO NGROUP_OCC_NEW = NGROUP_OCC_NEW + 1 NGROUP_REF(NGROUP_OCC_NEW) = N_CLASSES(IGRP1) NALTS = NALTS_OCCUP(IGRP) NALTS_OCCUP(NGROUP_OCC_NEW) = NALTS IF(NALTS.GT.0) THEN DO IE=1,NALTS NGROUP_ALT(IE,NGROUP_OCC_NEW)=N_CLASSES(NGROUP_ALT(IE,IGRP)) ENDDO ENDIF 400 CONTINUE ENDDO NGROUP_OCC = NGROUP_OCC_NEW cd STOP C C---Now find equavalence classes using relations RETURN END C SUBROUTINE OCCUPANCY_REFINE IMPLICIT NONE C C---Compile derivatives of function wrt occupancy groups. It works if C---all atoms divided into occupancy classe with their alternatives C C---It should be in the occup_flags.fh file INCLUDE 'atom_com.fh' INCLUDE 'pls_incl.fh' INCLUDE 'vitals.fh' INCLUDE 'occupancy_params.fh' C REAL AM,V,DV,WORKSPACE COMMON /LSQMTX/ AM( QQM),V( QQV),DV( QQV) REAL GX(3*MAXATOM),GQ(MAXATOM),GU1(6*MAXATOM) COMMON /REF_SPG/ GX,GU1,GQ COMMON /R_SCRATCH/ WORKSPACE(QQDEN) INTEGER IA,IQ,IGRP,IGRP1,IGRP2,NALTS,LV,NWORKSPACE,I,I1,I2, & NALTS1 C REAL DFDC(MAXALTS,MAXALT_GROUPS),SHIFTS_OCCUP(MAXALTS), & D2FDC2(MAXALTS,MAXALTS,MAXALT_GROUPS), & SHIFTS_GROUP(MAXALT_GROUPS*MAXALTS),SHIFTS_ALT(MAXALTS) LOGICAL ATOM_SELECT_OCCUP EXTERNAL ATOM_SELECT_OCCUP REAL TOLER,SUMSHIFTS DATA TOLER/1.0E-8/ C C---Initialise IQ = 0 WRITE(*,*)NGROUP_OCC,MAXALTS,MAXALT_GROUPS NWORKSPACE = QQDEN DO IGRP=1,NGROUP_OCC DO IGRP1=1,MAXALTS DFDC(IGRP1,IGRP) = 0.0 DO IGRP2=1,MAXALTS D2FDC2(IGRP1,IGRP2,IGRP) = 0.0 ENDDO ENDDO ENDDO C C---Now convert individual atomic occupancy derivatives to C---group occupancy derivatives DO IA=1,N_ATOM IF(ATOM_REF_FLAG(IA).GT.0) THEN IQ = IQ + 1 IF(IGRP_OCCUP(IA).GT.0) THEN IGRP1 = IGRP_OCCUP(IA) DO IGRP=1,NGROUP_OCC IF(IGRP1.EQ.NGROUP_REF(IGRP)) THEN C C---Group number coincideswith original group DFDC(1,IGRP) = DFDC(1,IGRP) + GQ(IQ) GOTO 100 ENDIF NALTS = NALTS_OCCUP(IGRP) IF(NALTS.GT.0) THEN DO IGRP2=1,NALTS C C--One of alternative groups IF(IGRP1.EQ.NGROUP_ALT(IGRP2,IGRP)) THEN IF(IGRP2.EQ.NALTS) THEN C C---Last group. It is not refined but used for constraints DO I=1,NALTS DFDC(I,IGRP) = DFDC(I,IGRP) - GQ(IQ) ENDDO ELSE C C---Other groups are refined as independent variables I1 = IGRP2+1 DFDC(I1,IGRP) = DFDC(I1,IGRP)+GQ(IQ) ENDIF GOTO 100 ENDIF ENDDO ENDIF ENDDO ENDIF ENDIF 100 CONTINUE ENDDO C C---Add lower triangle also DO IGRP=1,NGROUP_OCC NALTS = NALTS_OCCUP(IGRP) IF(NALTS.GT.1) THEN DO I1=1,NALTS_OCCUP(IGRP)-1 DO I2=I1+1,NALTS_OCCUP(IGRP) D2FDC2(I2,I1,IGRP) = D2FDC2(I1,I2,IGRP) ENDDO ENDDO ENDIF ENDDO C C----Solve equations with eigen value filtering. It is C----in principle unnecessary. Especially of number of alternatives C----is 2. Then only one parameter is refinable. DO IGRP=1,NGROUP_OCC NALTS = NALTS_OCCUP(IGRP) NALTS1 = NALTS IF(NALTS.LE.0) NALTS1 = 1 CALL EIGEN_FILTER_R(TOLER,D2FDC2(1,1,IGRP),NALTS1,MAXALTS, & DFDC(1,IGRP),SHIFTS_ALT,WORKSPACE,NWORKSPACE) I1 = NGROUP_REF(IGRP) SHIFTS_GROUP(I1) = SHIFTS_ALT(1) SUMSHIFTS = SHIFTS_ALT(1) IF(NALTS.GT.0) THEN DO I=1,NALTS-1 SHIFTS_GROUP(NGROUP_ALT(I,IGRP)) = SHIFTS_ALT(I+1) SUMSHIFTS = SUMSHIFTS + SHIFTS_ALT(I+1) ENDDO SHIFTS_GROUP(NGROUP_ALT(NALTS,IGRP)) = -SUMSHIFTS ENDIF ENDDO cd STOP LV = NVPOS + NVTMP C C---Now convert shifts of groups to those of atomic occupancies IQ = 0 DO IA=1,N_ATOM IF(ATOM_REF_FLAG(IA).GT.0) THEN IQ = IQ + 1 LV = LV + 1 DV(LV) = 0.0 IF(IGRP_OCCUP(IA).GT.0) THEN DV(LV) = SHIFTS_GROUP(IGRP_OCCUP(IA)) ELSE C C---Other atoms IF(ATOM_SELECT_OCCUP(IA)) THEN ENDIF ENDIF ENDIF ENDDO C RETURN END C LOGICAL FUNCTION ATOM_SELECT_OCCUP(IA) IMPLICIT NONE INCLUDE 'atom_com.fh' INTEGER IA ATOM_SELECT_OCCUP = .FALSE. IF(RES_NAME(I_RESID(IA))(1:3).EQ.'HOH') THEN ATOM_SELECT_OCCUP = .TRUE. ENDIF RETURN END C SUBROUTINE EQUIVALENCE_CLASSES(N_EL,N_EQUIV,N_RELATE1,N_RELATE2, & N_CLASSES) IMPLICIT NONE C C---Find class numbers using pairwise equivalence relations. First build C---tree then find class numbers using that tree. INTEGER N_EL,N_EQUIV INTEGER N_RELATE1(*),N_RELATE2(*),N_CLASSES(*) C INTEGER I,J,K INTEGER N_CLASSES1(1000) C C---Initialise DO I=1,N_EL N_CLASSES(I) = I ENDDO C DO I=1,N_EQUIV J = N_RELATE1(I) 1 CONTINUE IF(N_CLASSES(J).NE.J) THEN J = N_CLASSES(J) GOTO 1 ENDIF K = N_RELATE2(I) 2 CONTINUE IF(N_CLASSES(K).NE.K) THEN K = N_CLASSES(K) GOTO 2 ENDIF C C---This step is important for future searches. IF(J.GT.K) THEN N_CLASSES(J) = K ELSEIF(J.LT.K) THEN N_CLASSES(K) = J ENDIF ENDDO C C---Now walk thrhough upwards and find root of all trees. DO I=1,N_EL 3 CONTINUE IF(N_CLASSES(I).NE.N_CLASSES(N_CLASSES(I))) THEN N_CLASSES(I) = N_CLASSES(N_CLASSES(I)) GOTO 3 ENDIF ENDDO RETURN END C SUBROUTINE TREE_2_CLASS_NUMBER(N_EL,N_CLASSES,N_CLASS_N) IMPLICIT NONE INTEGER N_EL INTEGER N_CLASSES(*),N_CLASS_N(*) C INTEGER NUMBER_C,I,J,K C C---Find the class numbers in an increasing order. NUMBER_C = 1 N_CLASS_N(1) = NUMBER_C DO I=2,N_EL DO J=1,I-1 IF(N_CLASSES(I).EQ.N_CLASSES(J)) THEN N_CLASS_N(I) = N_CLASS_N(J) GOTO 10 ENDIF ENDDO NUMBER_C = NUMBER_C + 1 N_CLASS_N(I) = NUMBER_C 10 CONTINUE ENDDO cd WRITE(*,*)NUMBER_C,N_EL C RETURN END C subroutine d2da_radial(nref,nind,d2df,sigo) implicit none include 'weights.fh' integer nref integer nind(*) real d2df(*),sigo(*) if(sigma_refine_style.eq.'BINS') then call d2da_radial_bin(nref,nind,d2df,sigo) else call d2da_radial_exps(nref,nind,d2df,sigo) endif return end C SUBROUTINE D2DA_RADIAL_BIN(NREF,NIND,D2DF,SIGO) C IMPLICIT NONE C C---this subroutien find average values of D2DF in resolution bins. C---Some more things like epsilon symmetry, centricity should C---also be added INTEGER NREF INTEGER NIND(*) REAL D2DF(*),SIGO(*) INCLUDE 'weights.fh' INCLUDE 'celsym.fh' INCLUDE 'rharvest.fh' C INTEGER IR,IHH(3),ISYSAB,ICENT REAL EPSI,RSQ,RHO,STL_C REAL RSQ4,RSQ42,RSQ43,RSQ44,rsq3 REAL*8 WT1 C C--Test purposes INTEGER I_INTER REAL S0,DS C c---External functions and subprutines INTEGER IC INTEGER NWORKSPACE REAL*8 TOLER REAL*8 WORKSPACE(10000) REAL LSTLSQ,D2DF_WEIGHT C INTEGER IBIN,IBIN_S,NBIN_RAD1 INTEGER I,J,K INTEGER MAXBIN1 INTEGER NCYCLE_TCH,ITRY,N_TRY INTEGER N_USED_HERE REAL*8 DELTA,FVALUE,FVALUE_OLD REAL*8 DD1,DD2,D111 REAL*8 SHIFT_MAX,DEXP_D2 REAL DFDSIGMA,D2FDSIGMA2 REAL SEC_DER_C REAL*8 DFDS(MAXBIN+1),D2FDS2(MAXBIN+1,MAXBIN+1) REAL*8 SHIFTS(MAXBIN+1) real SEC_DER_DERS(MAXBIN+1) cd REAL*8 D_INV(40,40) real diff REAL TEMP1,TEMP2 cd REAL SCALE_SEC,B_SEC C EXTERNAL LSTLSQ,UNPACK,D2DF_WEIGHT C NWORKSPACE = 10000 MAXBIN1 = MAXBIN + 1 NBIN_ML1 = NBIN_ML + 1 NBIN_RAD1 = NBIN_RAD + 1 TOLER = 0.5E-8 cd WRITE(*,*) cd IF(NCYCLE_OVERALL.LE.1) THEN DO IBIN=1,NBIN_RAD+1 SEC_DER_BIN(IBIN) = 0.0 NREF_SEC_BIN(IBIN) = 0 ENDDO IBIN_S = 1 DO IR = 1,NREF IF(SIGO(IR).GT.0.0) THEN CALL INDTORS(NIND(IR),RSQ) CALL UNPACK(NIND(IR),IHH(1),IHH(2),IHH(3)) CALL CENTR(IHH,ICENT) CALL EPSLON(IHH,EPSI,ISYSAB) STL_C = SQRT(RSQ)**3 IF(.NOT.(STL_C.GE.SMEANB_RAD(IBIN_S).AND. & STL_C.LE.SMEANB_RAD(IBIN_S+1))) THEN DO IBIN=1,NBIN_RAD IF(STL_C.GT.SMEANB_RAD(IBIN).AND. & STL_C.LE.SMEANB_RAD(IBIN+1)) & THEN IBIN_S = IBIN GOTO 50 ENDIF ENDDO GOTO 60 ENDIF 50 CONTINUE DD1 = D2DF(IR)*EPSI*FLOAT(1+ICENT)/FLOAT(NSMULT) cd DD1 = 2.0*D2DF(IR)**2*EPSI*FLOAT(ICENT+1) cd DD1 = D2DF(IR) cd*EPSI*FLOAT(ICENT+1) SEC_DER_BIN(IBIN_S) = SEC_DER_BIN(IBIN_S) + DD1 NREF_SEC_BIN(IBIN_S) = NREF_SEC_BIN(IBIN_S) + 1 ENDIF 60 CONTINUE ENDDO DO IBIN=1,NBIN_RAD SEC_DER_BIN(IBIN) = alog(SEC_DER_BIN(IBIN)/NREF_SEC_BIN(IBIN)) cd WRITE(*,*)SEC_DER_BIN(IBIN),NREF_SEC_BIN(IBIN) ENDDO cd return C C--- Start refining parameters of smoothining function C TEMP1 = SEC_DER_BIN(1) SEC_DER_BIN(NBIN_RAD1) = SEC_DER_BIN(NBIN_RAD) DO IBIN=2,NBIN_ML TEMP2 = SEC_DER_BIN(IBIN) SEC_DER_BIN(IBIN) = (TEMP1+TEMP2)/2.0 TEMP1 = TEMP2 ENDDO cd SEC_DER_BIN(NBIN_ML1) = TEMP2 cd return DO IBIN=1,NBIN_RAD cd WRITE(*,*)'sec_der',smeanb_rad(ibin), cd & SMEANB_ML(IBIN),SEC_DER_BIN(IBIN) cd SEC_DER_BIN(IBIN) = ALOG(SEC_DER_BIN(IBIN)) ENDDO cd ENDIF c stop cd cd goto 100 if(nbin_rad.eq.1) goto 100 cd DO IC=1,5 DO I=1,NBIN_RAD1 DFDS(I) = 0.0 DO J=1,NBIN_RAD1 D2FDS2(I,J) = 0.0 ENDDO ENDDO C fvalue = 0.0 DO IR=1,NREF IF(SIGO(IR).GT.0.0) THEN CALL INDTORS(NIND(IR),RSQ) CALL UNPACK(NIND(IR),IHH(1),IHH(2),IHH(3)) CALL CENTR(IHH,ICENT) CALL EPSLON(IHH,EPSI,ISYSAB) rsq3 = sqrt(rsq**3) CALL SMOOTH_GAUSS_D(NBIN_RAD1,KERNEL_G_RAD,SMEANB_RAD, & SEC_DER_BIN,rsq3,SEC_DER_C,SEC_DER_DERS) cd write(10,*)rsq3,(smeanb_rad(i),sec_der_ders(i),i=1,nbin_ml1) dd2 = DBLE(AMAX1(-60.0,AMIN1(120.0,SEC_DER_C))) cd DD1 = 2.0*(D2DF(IR)**2*EPSI*FLOAT(ICENT+1)) cd DD1 = 2.0*D2DF(IR)*EPSI*FLOAT(ICENT+1) DD1 = DLOG(DBLE(D2DF(IR)*EPSI*FLOAT(ICENT+1)/FLOAT(NSMULT))) diff = (dd2 - dd1) fvalue = fvalue + diff**2 cd DFDSIGMA = -(DD1/SEC_DER_C-1.0/FLOAT(1+ICENT)) cd D2FDSIGMA2 = DFDSIGMA*DFDSIGMA cd D2FDSIGMA2 = DD1/SEC_DER_C cd D2FDSIGMA2 = DD1/SEC_DER_C DO I=1,NBIN_RAD1 DFDS(I) = DFDS(I) + diff*SEC_DER_DERS(I) D111 = SEC_DER_DERS(I) DO J=1,NBIN_RAD1 D2FDS2(I,J) = D2FDS2(I,J) + d111*SEC_DER_DERS(J) ENDDO ENDDO ENDIF ENDDO cd write(*,*)' fvalue ',fvalue C cd do i=1,nbin_ml1 cd write(*,*)dfds(i),(D2FDS2(I,j),j=1,nbin_ml1) cd enddo CALL DEIGEN_FILTER_R(TOLER,D2FDS2,NBIN_RAD1,MAXBIN1, & DFDS,SHIFTS,WORKSPACE,NWORKSPACE) DO IBIN=1,NBIN_RAD1 SEC_DER_BIN(IBIN) = SEC_DER_BIN(IBIN) - SHIFTS(IBIN) cd write(*,*) sec_der_bin(ibin) ENDDO cd write(*,*) cd ENDDO 100 continue C C---Tabulate for easy access C s0 = SMEANB_RAD(1) C cd OPEN(13,file='13.dat') cd write(*,*)'File opened ',N_POINTS_TABLE_S DO I=1,N_POINTS_TABLE_S+1 CALL SMOOTH_GAUSS(NBIN_RAD1,KERNEL_G_RAD,SMEANB_RAD,SEC_DER_BIN, & s0,SEC_DER_C) SEC_TABLE_SMOOTH(I) = dexp(dble(SEC_DER_C)) cd WRITE(13,*)s0,SEC_DER_C s0 = s0 + DELTA_S_TABLE ENDDO return cd stop write(*,*) DO I=1,NBIN_RAD1 WRITE(*,*)SEC_DER_BIN(I) ENDDO OPEN(13,file='13.dat') DO IR=1,NREF IF(SIGO(IR).GT.0.0) THEN CALL INDTORS(NIND(IR),RSQ) CALL UNPACK(NIND(IR),IHH(1),IHH(2),IHH(3)) CALL CENTR(IHH,ICENT) CALL EPSLON(IHH,EPSI,ISYSAB) CALL SMOOTH_GAUSS(NBIN_RAD1,KERNEL_G_RAD,SMEANB_RAD, & SEC_DER_BIN,sqrt(RSQ**3),SEC_DER_C) DD1 = D2DF(IR)*EPSI*FLOAT(1+ICENT)/FLOAT(NSMULT) cd DD1 = 2.0*(D2DF(IR)**2*EPSI*FLOAT(ICENT+1)) WRITE(13,*)RSQ,DD1,exp(SEC_DER_C) ENDIF ENDDO C C--- Tabulate now for future use cd return stop C RETURN END C SUBROUTINE D2DA_RADIAL_EXPS(NREF,NIND,D2DF,SIGO) C IMPLICIT NONE C C---this subroutien find average values of D2DF in resolution bins. C---Some more things like epsilon symmetry, centricity should C---also be added INTEGER NREF INTEGER NIND(*) REAL D2DF(*),SIGO(*) INCLUDE 'weights.fh' INCLUDE 'celsym.fh' INCLUDE 'rharvest.fh' C INTEGER IR,IHH(3),ISYSAB,ICENT REAL EPSI,RSQ,RHO,STL_C REAL RSQ4,RSQ42,RSQ43,RSQ44,rsq3 REAL*8 WT1 C C--Test purposes INTEGER I_INTER REAL S0,DS C c---External functions and subprutines INTEGER IC INTEGER NWORKSPACE REAL*8 TOLER/1.0D-7/ REAL*8 WORKSPACE(300) REAL LSTLSQ,D2DF_WEIGHT C INTEGER IBIN,IBIN_S INTEGER I,J,K INTEGER NCYCLE_TCH,ITRY,N_TRY INTEGER N_USED_HERE REAL*8 DELTA,FVALUE,FVALUE_OLD REAL DD1,DD2 REAL SHIFT_MAX,DEXP_D2 REAL SEC_DER_C integer MAXCOEF PARAMETER (MAXCOEF = 6) REAL*8 DFDS(MAXCOEF),D2FDS2(MAXCOEF,MAXCOEF) REAL*8 SHIFTS(MAXCOEF) cd REAL SCALE_SEC,B_SEC C EXTERNAL LSTLSQ,UNPACK,D2DF_WEIGHT C NWORKSPACE = 100 DO I=1,4 DFDS(I) = 0.0 DO J=1,4 D2FDS2(I,J) = 0.0 ENDDO ENDDO C DO IR=1,NREF IF(SIGO(IR).GT.0.0) THEN CALL INDTORS(NIND(IR),RSQ) CALL UNPACK(NIND(IR),IHH(1),IHH(2),IHH(3)) CALL CENTR(IHH,ICENT) CALL EPSLON(IHH,EPSI,ISYSAB) rsq4 = rsq/4.0 rsq42 = rsq4*rsq4 rsq43 = rsq42*rsq4 cd dd2 = scale_sec-b_sec*rsq4+b1_sec*rsq42+b2_sec*rsq43 DD1 = alog(D2DF(IR)*EPSI*FLOAT(ICENT+1)/FLOAT(NSMULT)) dfds(1) = dfds(1) + dd1 dfds(2) = dfds(2) - dd1*rsq4 dfds(3) = dfds(3) + dd1*rsq42 dfds(4) = dfds(4) + dd1*rsq43 d2fds2(1,1) = d2fds2(1,1) + 1.0 d2fds2(1,2) = d2fds2(1,2) - rsq4 d2fds2(2,2) = d2fds2(2,2) + rsq4*rsq4 d2fds2(1,3) = d2fds2(1,3) + rsq42 d2fds2(2,3) = d2fds2(2,3) - rsq42*rsq4 d2fds2(3,3) = d2fds2(3,3) + rsq42*rsq42 d2fds2(1,4) = d2fds2(1,4) + rsq43 d2fds2(2,4) = d2fds2(2,4) - rsq43*rsq4 d2fds2(3,4) = d2fds2(3,4) + rsq43*rsq42 d2fds2(4,4) = d2fds2(4,4) + rsq43*rsq43 ENDIF ENDDO d2fds2(2,1) = d2fds2(1,2) d2fds2(3,1) = d2fds2(1,3) d2fds2(3,2) = d2fds2(2,3) d2fds2(4,1) = d2fds2(1,4) d2fds2(4,2) = d2fds2(2,4) d2fds2(4,3) = d2fds2(3,4) CALL DEIGEN_FILTER_R(TOLER,D2FDS2,4,MAXCOEF, & DFDS,SHIFTS,WORKSPACE,NWORKSPACE) scale_sec = shifts(1) b_sec = shifts(2) b1_sec = shifts(3) b2_sec = shifts(4) 100 continue C C---Tabulate for easy access C s0 = SMEANB_RAD(1) C cd OPEN(13,file='13.dat') cd write(*,*)'File opened ',N_POINTS_TABLE_S DO I=1,N_POINTS_TABLE_S+1 rsq4 = s0**(2.0/3.0)/4.0 sec_table_smooth(i) = exp(scale_sec-rsq4*(b_sec- & rsq4*(b1_sec+rsq4*b2_sec))) s0 = s0 + DELTA_S_TABLE ENDDO return cd stop cd write(*,*) DO I=1,NBIN_ML1 WRITE(*,*)SEC_DER_BIN(I) ENDDO OPEN(13,file='13.dat') i_inter = -1 DO IR=1,NREF IF(SIGO(IR).GT.0.0) THEN CALL INDTORS(NIND(IR),RSQ) CALL UNPACK(NIND(IR),IHH(1),IHH(2),IHH(3)) CALL CENTR(IHH,ICENT) CALL EPSLON(IHH,EPSI,ISYSAB) s0 = sqrt(rsq) sec_der_c = d2df_weight(s0,i_inter) DD1 = D2DF(IR)*EPSI*FLOAT(1+ICENT)/FLOAT(NSMULT) WRITE(13,*)RSQ,DD1,exp(SEC_DER_C) ENDIF ENDDO C C--- Tabulate now for future use cd return stop C RETURN END C REAL*8 FUNCTION D2DF_WEIGHT(S_CURRENT,I_INTER) IMPLICIT NONE C c--Calculates current value of second derivative at the given point (radial part) C--This version uses linear interpolation. It will have to be checked and probably C--other interpolation techniques will have to be used C INTEGER I_INTER REAL H_SM REAL S_CURRENT INCLUDE 'weights.fh' INCLUDE 'celsym.fh' INTEGER I_ACCESS REAL*8 CURRENT_VALUE REAL CURRENT_VALUE1,CURRENT_VALUE2 REAL S21,S22,S23 REAL S_CURRENT2 REAL SIGMA_IN1,SIGMA_IN,SCURRENT EXTERNAL LINTER_VALUE2,GAUSS_SMOOTH C cd IF(I_INTER.LE.0) THEN cd H_SM = -0.1 cd I_INTER = 1 cd ENDIF cd CALL GAUSS_SMOOTH(H_SM,NBIN_RAD,SMEANB_RAD,D2DF_ML_RAD_BIN, cd & S_CURRENT,CURRENT_VALUE) cd IF(S_RAD.LE.0.0.AND.ABS(B_RAD).GT.2.0) THEN cd CALL LINTER_VALUE2(NBIN_ML,SMAXB_ML,SEC_DER_BIN, cd & S_CURRENT,I_INTER,CURRENT_VALUE) cd CURRENT_VALUE = CUR1 cd SCURRENT = S_CURRENT/SMAX_TCH cd CALL CALC_TCHEB_APPROX_CLENSHW(N_COEFS,SCURRENT,TCH_COEFS, cd & SIGMA_IN) cd CALL CALC_TCHEB_APPROX_CLENSHW(N_COEFS,SCURRENT, cd & TCH_COEFS_SIGMA1,SIGMA_IN1) cd ELSE cd CURRENT_VALUE = S_RAD*EXP(-B_RAD*S_CURRENT**2/4.0) cd ENDIF IF(SIGMA_REFINE_STYLE.EQ.'BINS') THEN S_CURRENT2 = S_CURRENT**3 I_ACCESS = NINT((S_CURRENT2-SMEANB_RAD(1))/DELTA_S_TABLE)+1 I_ACCESS = MAX(1,MIN(N_POINTS_TABLE_S+1,I_ACCESS)) CURRENT_VALUE2 = SEC_TABLE_SMOOTH(I_ACCESS) D2DF_WEIGHT = DBLE(CURRENT_VALUE2) ELSE S21 = (S_CURRENT/2.0)**2 S22 = S21*S21 S23 = S22*S21 D2DF_WEIGHT = DEXP(DBLE(SCALE_SEC-B_SEC*S21+B1_SEC*S22+ & B2_SEC*S23)) ENDIF cd CALL LINTER_VALUE2(NBIN_ML,SMEANB_ML,SEC_DER_BIN, cd & S_CURRENT2,I_INTER,CURRENT_VALUE2) cd CALL SMOOTH_GAUSS(NBIN_ML1,KERNEL_G,SMEANB_ML,SEC_DER_BIN, cd & S_CURRENT2,CURRENT_VALUE2) cd S21 = (S_CURRENT/2.0)**2 cd S22 = S21*S21 cd S23 = S22*S21 cd CURRENT_VALUE = DEXP(DBLE(SCALE_SEC-B_SEC*S21+B1_SEC*S22+ cd & B2_SEC*S23)) cd & 2.0*(SIGMA_IN+SIGMA_ML_SCALE_OVER- cd & SIGMA_ML_B_OVER*S_CURRENT**2/4.0))) cd WRITE(*,*)D2DF_WEIGHT,CURRENT_VALUE2,S_CURRENT2,SMEANB_ML(1), cd & SMEANB_ML(NBIN_ML1) cd/DBLE(NSMULT) END C subroutine write_reflections_r(nref,fo,sigo,fc,phase,nind) C implicit none include 'refi_flags.fh' integer nref integer nind(nref) real fo(*),sigo(*) real fc(*),phase(*) C integer iout_unit,ifail,ll integer npart1 integer ir,ip integer hkl(3) real ss character file_name*512 real lstlsq external lstlsq,unpack C call ugtenv('SFOUT',file_name) if(file_name(1:1).eq.' ') return npart1 = npart + 1 iout_unit = 0 ifail = -1 ll = 0 call ccpdpn(iout_unit,'SFOUT','UNKNOWN','F',ll,ifail) do ir=1,nref call unpack(nind(ir),hkl(1),hkl(2),hkl(3)) ss = lstlsq(1,hkl(1),hkl(2),hkl(3)) write(iout_unit,*)hkl,ss,(fc(ir+nobs*ip), & phase(ir+nobs*ip),ip=0,npart) enddo stop return end