diff -u -N -r xplor-nih-2.9.2/source/Makefile PARAxplor-nih-2.9.2/source/Makefile --- xplor-nih-2.9.2/source/Makefile 2003-02-13 17:00:24.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/Makefile 2003-12-05 14:34:36.000000000 +0100 @@ -83,7 +83,9 @@ xrama.f trio_faster_newderiv.f susc_anis.f \ rama_gaussians.f orientations.f xdipolar_coup.f \ diff_anis.f anisotropy.f angledb.f \ - vectangl.f + vectangl.f \ + fantacross.f fantalin.f xangle.f xccr.f xfantaccr.f \ + xfantadipo.f xdipo_pcs.f xdipo_rdc.f xt1dist.f SERIAL_SRC += hb_dist_angle.f diff -u -N -r xplor-nih-2.9.2/source/ener.fcm PARAxplor-nih-2.9.2/source/ener.fcm --- xplor-nih-2.9.2/source/ener.fcm 2003-01-06 20:27:05.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/ener.fcm 2003-12-05 14:34:36.000000000 +0100 @@ -2,7 +2,12 @@ C C total number of energy terms (note the free slots!) INTEGER NENERT - PARAMETER (NENERT=52) +C ------------------------------------------------------------------ +C Number of energy terms updated to deal with PARArestraints (2003) +C +C PARAMETER (NENERT=52) + PARAMETER (NENERT=56) +C ------------------------------------------------------------------ C C DSWGHT weight of the individual double selections in the Hamiltonian. C DSAVWT weight of the individual double selections in the @@ -26,6 +31,11 @@ INTEGER SSDIPO, SSDANI, SSXDIP, SSXRAM, SSVEAN, SSDCSA INTEGER SSCRIPT, SSHBDA, SSPMAG INTEGER SSGBSELF, SSGBINT, SSPGBSLF, SSPGBINT +C ------------------------------------------------------------------ +C These are the energy terms related to PARArestraints (2003) +C + INTEGER SSXPCS, SSXRDC, SSXCCR, SSXANGLE +C ------------------------------------------------------------------ C C C @@ -44,6 +54,11 @@ PARAMETER (SSCRIPT=46, SSHBDA=47) PARAMETER (SSGBSELF= 48, SSPGBSLF=49, SSGBINT=50, SSPGBINT=51) PARAMETER (SSPMAG= 52) +C ------------------------------------------------------------------ +C These are the energy terms related to PARArestraints (2003) +C + PARAMETER (SSXPCS=53, SSXRDC=54, SSXCCR=55, SSXANGLE=56) +C ------------------------------------------------------------------ C C >>>>>>> NOTE: the GB self terms MUST be numbered consecutively <<<<<<<<< C (because they are calculated first and must not be reinitialized by ENERG2) diff -u -N -r xplor-nih-2.9.2/source/energy.f PARAxplor-nih-2.9.2/source/energy.f --- xplor-nih-2.9.2/source/energy.f 2003-09-30 23:48:05.000000000 +0200 +++ PARAxplor-nih-2.9.2/source/energy.f 2003-12-05 14:34:36.000000000 +0100 @@ -8,6 +8,14 @@ include 'ener.fcm' include 'deriv.fcm' include 'timer.fcm' + +c ------------------------------------------------------------------ +c FANTALIN stuff (2003) +c + double precision fenergtot + common /fenerg/fenergtot +c ------------------------------------------------------------------ + C parameters DOUBLE PRECISION ZERO PARAMETER (ZERO=0.0D0) @@ -46,6 +54,12 @@ RENR(SSENER)=TOTAL CALL DECLAR( ANER(SSENER), 'DP', ' ', DCVAL, RENR(SSENER) ) C +c ------------------------------------------------------------------ +c FANTALIN stuff (2003) +c + fenergtot=total +c ------------------------------------------------------------------ + C compute gradient and put in the accumulation slot SSSD TOTAL=ZERO NDIM=0 @@ -224,6 +238,34 @@ RENR(SSDIPO) = RENR(SSXDIP) IF (TIMER.GT.1) CALL DSPCPU(' ENERGY: after Dipolar term') ENDIF +c ------------------------------------------------------------------ +c Energy terms related to PARArestraints (2003) +c +C +C pseudocontact shifts energy + IF (QENER(SSXPCS)) THEN + CALL EXPCS(RENR(SSXPCS),'ENERGY ') + IF (TIMER.GT.1) CALL DSPCPU(' ENERGY: after PCS term') + ENDIF +C +C residual dipolar couplings energy + IF (QENER(SSXRDC)) THEN + CALL EXRDC(RENR(SSXRDC),'ENERGY ') + IF (TIMER.GT.1) CALL DSPCPU(' ENERGY: after RDC term') + ENDIF +C +C cross-correlation rates energy + IF (QENER(SSXCCR)) THEN + CALL EXCCR(RENR(SSXCCR),'ENERGY ') + IF (TIMER.GT.1) CALL DSPCPU(' ENERGY: after CCR term') + ENDIF +C +C RDCANGLES energy + IF (QENER(SSXANGLE)) THEN + CALL EXANGLES(RENR(SSXANGLE),'ENERGY ') + IF (TIMER.GT.1) CALL DSPCPU(' ENERGY: after ANG term') + ENDIF +c ------------------------------------------------------------------ C C Susceptibility anisotropy energy IF(QENER(SSSANI)) THEN @@ -440,6 +482,14 @@ ANER(SSDANI)='DANI' ANER(SSXDIP)='XDIP' ANER(SSDIPO)='DIPO' +c ------------------------------------------------------------------ +c Slots for PARArestraints (2003) +c + ANER(SSXPCS)='XPCS' + ANER(SSXRDC)='XRDC' + ANER(SSXCCR)='XCCR' + ANER(SSXANGLE)='XANG' +c ------------------------------------------------------------------ ANER(SSDCSA)='DCSA' ANER(SSHBDA)='HBDA' ANER(SSVEAN)='VEAN' @@ -518,6 +568,11 @@ &' HBON|XREF|NOE|PVDW|PELE|HARM|CDIH|USER|SBOU|NCS|PLAN|', &' PGBS|PGBI|COUP|CARB|PROT|ONEJ|RAMA|XRAM|ANDB|ANIS|TRIO', &' ORIE|SANI|DANI|DIPO|DCSA|HBDA|VEAN|SCRIpt', +c ------------------------------------------------------------------ +c Energy terms related to PARArestraints (2003) +c + &' XPCS|XRDC|XCCR|XANG', +c ------------------------------------------------------------------ &' remark: wildcards are allowed for ' C-DOCUMENTATION-SOURCE-END C diff -u -N -r xplor-nih-2.9.2/source/fantacross.f PARAxplor-nih-2.9.2/source/fantacross.f --- xplor-nih-2.9.2/source/fantacross.f 1970-01-01 01:00:00.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/fantacross.f 2003-12-05 14:34:36.000000000 +0100 @@ -0,0 +1,482 @@ +C ------------------------------------------------------------------ + SUBROUTINE FANTALINCCR(DJI,DJK,ANGFACT,NNAT,XOBS,XTOLPROT, + & NOMEFILE1,ERWPROT,ERRORE) +C ------------------------------------------------------------------ +C Handles the stuff for the calculation of CCR constant coefficient +C +C CX,CY,CZ(STRUCTURE,ATOMNUMBER) keep coordinates of nuclei +C FX,FY,FZ(STRUCTURE,ATOMNUMBER) keep coordinates of metal +C NAT is a vector with the number of atoms +C OBS is a vector with the observed cross correlation rates +C TOLPROT is a vector with the errors on CCR +C NOMEFILE1 is the name of the file to save deviations +C ERWPROT is a vector with the weights +C ERRORE switches the use(1) of weights +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + INCLUDE 'supccr.fcm' + INCLUDE 'fantaxplor.fcm' + + INTEGER NNAT,METHOD + DOUBLE PRECISION AK,A + COMMON /CRUN/AK + DOUBLE PRECISION DJI(1000),DJK(1000),ANGFACT(1000), + & XOBS(1000),XTOLPROT(1000),ERWPROT(1000) + CHARACTER NOMEFILE1*132 + INTEGER ERRORE + + DO I=1,1000 + AK1(I)=ANGFACT(I) + RK1(I)=DJI(I) + RK2(I)=DJK(I) + END DO + + PRINT*,' ' + PRINT*,'===---FANTACROSS---===' +C ------------------------------------------------------------------ +C IHP is the number of CCR +C ------------------------------------------------------------------ + NAT=NNAT + IHP=NNAT + NHP=NNAT + PRINT*,'IHP is:',IHP +C ------------------------------------------------------------------ +C NFE is the number of paramagnetic centers +C ------------------------------------------------------------------ + NFE=1 + NIA=0 +C ------------------------------------------------------------------ +C Reference system not active +C ------------------------------------------------------------------ + NSYSTEM=0 +C ------------------------------------------------------------------ +C Calculation is done on a single structure +C ------------------------------------------------------------------ + NSTR=1 +C ------------------------------------------------------------------ +C Errors are set to 0 +C ------------------------------------------------------------------ + PERC=0 +C ------------------------------------------------------------------ +C Grid setting +C ------------------------------------------------------------------ + NGRID=1 +C ------------------------------------------------------------------ +C Weights and multiplicity are set to 1 +C ------------------------------------------------------------------ + DO NIA=1,NAT + MLPROT(NIA)=1 + IF(ERRORE.EQ.0) THEN + WPROT(NIA)=1 + ELSE + WPROT(NIA)=ERWPROT(NIA) + END IF + END DO + + OPEN(322,FILE='FANTACROSS.LOG') +C ------------------------------------------------------------------ +C Vectors are passed from Xplor-NIH core +C ------------------------------------------------------------------ + DO NIA=1,NAT + OBS(NIA)=XOBS(NIA) +C ------------------------------------------------------------------ +C No tolerance in fitting +C ------------------------------------------------------------------ + TOLPROT(NIA)=0 + END DO + CLOSE (322) + FILENAME3=NOMEFILE1 + PRINT*,'NOMEFILE1 is:',NOMEFILE1 + OPEN(744,FILE=FILENAME3) + + CALL CSIMPLEXRUN() + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE CSIMPLEXRUN() +C ------------------------------------------------------------------ +C Calculates the CCR constant coefficient which minimizes the +C experimental-to-calculated squared difference +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + INCLUDE 'supccr.fcm' + + DIMENSION SIMP(MP,NP),EXTR(NP),VAL(MP) + DIMENSION ZRR(3,3) + + OLDRESID=1.E+9 + RESID=2.E+9 + IVIOLATION=2000000 + IOLDVIO=1000000 +C ------------------------------------------------------------------ +C Initial conditions grid +C ------------------------------------------------------------------ + DELTA=3.14159/DBLE(NGRID+1) + DO 1000 L=1,NGRID + DO 999 NK=1,NFE + SIMP(1,(NK-1)*MPAR+1)=1000.*(1.-2*RAND()) +999 CONTINUE + DO 10 I=1,MP + DO 10 J=1,NP + SIMP(I,J)=SIMP(1,J) + IF (I.EQ.1) THEN + SIMP(1,J)=SIMP(1,J)*0.8 + ENDIF + IF (I.EQ.J+1) THEN + SIMP(I,J)=SIMP(I,J)/0.8*1.2 + ENDIF +10 CONTINUE + DO 11 I=1,MP + DO 12 J=1,NP +12 EXTR(J)=SIMP(I,J) + VAL(I)=CCASH(EXTR) +11 CONTINUE + CALL CSIMPCALC(SIMP,VAL) +C ------------------------------------------------------------------ +C Criterion of minimization of optimal solution +C ------------------------------------------------------------------ + IF (OLDRESID.GE.RESID) THEN + OLDRESID=RESID + IOLDVIO=IVIOLATION + DO 13 I=1,NFE + OPTPHI(I)=PHI(I) +13 CONTINUE + ENDIF +1000 CONTINUE + + DO 998 NK=1,NFE + EXTR((NK-1)*MPAR+1)=OPTPHI(NK) +998 CONTINUE + + PRINT*,' ' + WRITE(*,FMT='(1X,A)') ' ***** BEST SOLUTION *****' + PRINT*,' ' + WRITE(*,20) IOLDVIO,OLDRESID + WRITE(*,21) (OPTPHI(I),I=1,NFE) + WRITE(*,*) 'Value of K =',OPTPHI(1) + WRITE(744,*) 'Value of K =',OPTPHI(1) + +20 FORMAT (1X,' VIOLATIONS ',I6,1X,' RESIDUALS ',F10.3) +21 FORMAT (1X,'K = ',F15.3) + + ERR=CCASH(EXTR) + CALL CDISPSHIFT + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE CSIMPCALC(P,Y) +C ------------------------------------------------------------------ +C Simplex calculation +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + INCLUDE 'supccr.fcm' + INCLUDE 'fantaxplor.fcm' + + PARAMETER (AL=1.,BE=.5,GA=2.) + DIMENSION P(MP,NP),Y(MP),PR(MP),PT(MP),PE(MP), + & PC(MP),PK(MP),PBAR(MP),OLDP(NP),EXTR(MP) +C ------------------------------------------------------------------ +C Algorithm initialization +C ------------------------------------------------------------------ + ITER=0 + RNDOM=0 +1 PERM=1 + DO 1000 WHILE (PERM.GT.0) + PERM=0 + DO 1100 I=1,NP + IF (Y(I).GT.Y(I+1)) THEN + APP=Y(I+1) + Y(I+1)=Y(I) + Y(I)=APP + DO 1200 J=1,NP + APP=P(I+1,J) + P(I+1,J)=P(I,J) + P(I,J)=APP +1200 CONTINUE + PERM=PERM+1 + ENDIF +1100 CONTINUE +1000 CONTINUE + + AVERAGE=0 + ASTDEV=0 + DO 100 I=1,MP +100 AVERAGE=AVERAGE+Y(I) + AVERAGE=AVERAGE/MP + DO 101 I=1,MP +101 ASTDEV=ASTDEV+(Y(I)-AVERAGE)**2. + ERR=DSQRT(ASTDEV/MP) + RESID=Y(1) + IF (ITER.GT.39999.OR.ERR.LT.1.D-9) THEN + WRITE (*,200) ITER,IVIOLATION,RESID + WRITE (*,201) (P(1,(I-1)*MPAR+5),I=1,NFE) +200 FORMAT (1X,'Iters ',I6,1X,' Violations ',I6, + & 1X,'Residuals= ',F10.3) +201 FORMAT (1X,'K = ',F9.3) + DO 997 NK=1,NFE + PHI(NK)=P(1,(NK-1)*MPAR+1) +997 CONTINUE + RETURN + ENDIF + + ITER=ITER+1 + DO 12 J=1,NP +12 PBAR(J)=0 + DO 13 I=1,NP + DO 13 J=1,NP +13 PBAR(J)=PBAR(J)+P(I,J) + DO 14 J=1,NP + PBAR(J)=PBAR(J)/NP + PR(J)=(1.+AL)*PBAR(J)-AL*P(MP,J) +14 CONTINUE + YR=CCASH(PR) + DO 15 J=1,NP + PK(J)=PR(J) +15 CONTINUE + YK=YR + IF (YR.LT.Y(NP)) THEN + IF (YR.LT.Y(1)) THEN + DO 16 J=1,NP + PE(J)=GA*PR(J)+(1.-GA)*PBAR(J) +16 CONTINUE + YE=CCASH(PE) + ITER=ITER+1 + IF (YE.LT.Y(1)) THEN + DO 17 J=1,NP + PK(J)=PE(J) +17 CONTINUE + YK=YE + ENDIF + ENDIF + ELSE + DO 18 J=1,NP + PT(J)=P(MP,J) +18 CONTINUE + YT=Y(MP) + IF (YR.LT.YT) THEN + DO 19 J=1,NP + PT(J)=PR(J) +19 CONTINUE + YT=YR + ENDIF + DO 20 J=1,NP + PC(J)=BE*PT(J)+(1.-BE)*PBAR(J) +20 CONTINUE + YC=CCASH(PC) + ITER=ITER+1 + IF (YC.LT.Y(NP)) THEN + DO 21 J=1,NP + PK(J)=PC(J) +21 CONTINUE + YK=YC + ELSE + DO 22 I=1,NP + OLDP(I)=P(2,I) +22 CONTINUE + DO 24 I=2,NP + DO 23 J=1,NP + P(I,J)=.5*(P(1,J)+P(I,J)) +23 CONTINUE + PR(I)=P(I,I) + Y(I)=CCASH(PR) +24 CONTINUE + ITER=ITER+NP-1 + DO 25 I=1,NP + OLDP(I)=ABS(OLDP(I)-P(2,I)) +25 CONTINUE + OLDPOPPA=0 + DO 251 I=1,NP-1 + OLDPOPPA=OLDPOPPA+OLDP(I) +251 CONTINUE + IF (ABS(OLDPOPPA).LT.1.D-8.AND.RNDOM.LT.10) THEN + RNDOM=RNDOM+1 + WRITE(*,FMT='(BN,A,F16.6)') 'RANDOM ',Y(1) + DO 2001 IL=1,MP + DO 2002 IM=1,NP + RMR=1+RAND() + P(IL,IM)=P(IL,IM)*RMR +2002 CONTINUE +2001 CONTINUE + DO 2011 IL=1,MP + DO 2012 JM=1,NP +2012 EXTR(JM)=P(IL,JM) + Y(IL)=CCASH(EXTR) +2011 CONTINUE + GOTO 1 + ENDIF + IF (ABS(OLDPOPPA).LT.1.D-9) THEN + ITER=50000 + ENDIF + DO 26 J=1,NP + PK(J)=0.5*(P(1,J)+P(MP,J)) +26 CONTINUE + YK=CCASH(PK) + ITER=ITER+1 + ENDIF + ENDIF + DO 27 J=1,NP + P(MP,J)=PK(J) +27 CONTINUE + Y(MP)=YK + GOTO 1 + END +C ------------------------------------------------------------------ + FUNCTION CCASH(VETT) +C ------------------------------------------------------------------ +C Calculation of deviations and violations +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + INCLUDE 'supccr.fcm' + + DIMENSION VETT(NP) + INTEGER M1 + + IVIOLATION=0 + TMP1=0 + I=1 + DO WHILE (I.LE.NHP*NSTR) + CSHIFT(I)=0.0 + I=I+1 + ENDDO + + DO 1 N=1,NSTR + DO 2 M=1,NFE + P=VETT((M-1)*MPAR+1) + TMP2=0.0 + I=1 + IHP=(N-1)*NHP + DO 10 WHILE (I.LE.NHP) + NPROTML=MLPROT(IHP+I) + CSHIFT(IHP+I)=CSHIFT(IHP+I)+AK1(IHP+I)/(RK2(IHP+I)**3* + & RK1(IHP+I)**3) + I=I+1 +10 CONTINUE +2 CONTINUE + DO 3 I=1,NHP + TMP2=ABS(CSHIFT(IHP+I)-OBS(IHP+I))-TOLPROT(IHP+I) + IF (TMP2.GT.0.0) THEN + IVIOLATION=IVIOLATION+1 + TMP1=TMP1+TMP2**2*WPROT(IHP+I)/DBLE(MLPROT(IHP+I)) + ENDIF +3 CONTINUE +1 CONTINUE + + CCASH=TMP1 + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE CDISPSHIFT +C ------------------------------------------------------------------ +C Display of violations +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + INCLUDE 'supccr.fcm' + INCLUDE 'fantaxplor.fcm' + + CHARACTER FILESTR*21,SS*40 + DIMENSION VV(MAXSTR,MAXOS),VMAX(MAXOS),VMIN(MAXOS) + DIMENSION NNMAX(MAXOS),NNMIN(MAXOS) + + RM=0 + TMP1=0 + SMED=0 + STDEV=0 + DO 2 I=1,NSTR + DO 2 J=1,NHP + TMP2=ABS(OBS(J+(I-1)*NHP)-CSHIFT(J+(I-1)*NHP))- + & TOLPROT(J+(I-1)*NHP) + IF (TMP2.GT.0.0) THEN + IF (CSHIFT(J+(I-1)*NHP).GT.OBS(J+(I-1)*NHP)) THEN + VV(I,J)=TMP2 + ELSE + VV(I,J)=-TMP2 + ENDIF + ELSE + VV(I,J)=0. + ENDIF +2 CONTINUE + DO 44 I=1, NHP + VMAX(I)=0. + VMIN(I)=1000. +44 CONTINUE + DO 4 J=1,NHP + DO 4 I=1,NSTR + IF (ABS(VV(I,J)).GT.ABS(VMAX(J))) THEN + VMAX(J)=VV(I,J) + ENDIF + IF ((ABS(VV(I,J)).LT.ABS(VMIN(J))).AND.(VV(I,J).NE.0.)) THEN + VMIN(J)=VV(I,J) + ENDIF +4 CONTINUE + DO 7 J=1,NHP + DO 8 I=1,NSTR + IF (ABS(VV(I,J)).GT.0.) THEN + SS(I:I+1)="*" + ELSE + SS(I:I+1)=" " + ENDIF +8 CONTINUE + IF (VMIN(J).EQ.1000.) THEN + VMIN(J)=0. + ENDIF +7 CONTINUE + + WRITE(744,'(1X)') + WRITE(744,'(1X)') + WRITE(744,'(1X)') + WRITE(744,'(1X,A)') + & ' PROTONS OBS. CALC. ERR.^2' + + FILESTR='USELESS.FILE' + NIA=0 + REWIND(1) + DO 3 J=1,NHP*NSTR + TMP2=ABS(OBS(J)-CSHIFT(J))-TOLPROT(J) + NIA=NIA+1 + IF (TMP2.GT.0.0) THEN + IVIOLATION=IVIOLATION+1 + RM=RM+1 + SMED=SMED+TMP2**2*WPROT(J)/MLPROT(J) + TMP1=TMP2**2*WPROT(J)/MLPROT(J) + WRITE (744,'(1X,A4,1X,A4,1X,A4,1X,F7.3,1X,F7.3,1X,F10.5)') + & NRESIDFA(NIA),NRESFA(NIA),NTYPEFA(NIA), + & OBS(J),CSHIFT(J),TMP1 +15 FORMAT(1X,A,1X,F7.3,1X,F7.3,1X,F7.3) + ELSE + WRITE (744,'(1X,A4,1X,A4,1X,A4,1X,F7.3,1X,F7.3)') + & NRESIDFA(NIA),NRESFA(NIA),NTYPEFA(NIA), + & OBS(J),CSHIFT(J) + ENDIF +3 CONTINUE + SMED=SMED/RM + + DO 41 J=1,NHP*NSTR + TMP2=ABS(OBS(J)-CSHIFT(J))-TOLPROT(J) + IF (TMP2.GT.0.0) THEN + STDEV=STDEV+(TMP2**2-MED)**2*WPROT(J)/MLPROT(J) + ENDIF +41 CONTINUE + STDEV=STDEV/RM + + PRINT*,' ' + WRITE(744,'(A)') ' ' + WRITE(*,11) 'Mean value of error =',SMED +11 FORMAT(1X,A,1X,F10.6) + WRITE(744,12) 'Mean value of error =',SMED +12 FORMAT(1X,A,1X,F10.6) + WRITE(*,13) 'Standard deviation =',STDEV +13 FORMAT(1X,A,1X,F10.6) + WRITE(744,14) 'Standard deviation =',STDEV +14 FORMAT(1X,A,1X,F10.6) + CLOSE(744) + + RETURN + END diff -u -N -r xplor-nih-2.9.2/source/fantalin.f PARAxplor-nih-2.9.2/source/fantalin.f --- xplor-nih-2.9.2/source/fantalin.f 1970-01-01 01:00:00.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/fantalin.f 2003-12-05 14:34:36.000000000 +0100 @@ -0,0 +1,1680 @@ +C ------------------------------------------------------------------ + SUBROUTINE FANTALINEXE(XPX,XPY,XPZ,XPMX,XPMY,XPMZ,NNAT,XOBS, + & XTOLPROT,METHOD,NOMEFILE1,ERWPROT,ERRORE, + & LOCFLG) +C ------------------------------------------------------------------ +C Handles the stuff for the calculation of tensor parameters +C +C METHOD switches between PCS(1) and RDC(0) +C CX,CY,CZ(STRUCTURE,ATOMNUMBER) keep coordinates of nuclei +C FX,FY,FZ(STRUCTURE,ATOMNUMBER) keep coordinates of metal +C NAT is a vector with the number of atoms +C OBS is a vector with the observed pseudocontact shifts +C TOLPROT is a vector with the errors on PCS +C NOMEFILE1 is the name of the file to save deviations +C ERWPROT is a vector with the weights +C ERRORE switches the USE(1) of weights +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + INCLUDE 'sup2.fcm' + INCLUDE 'fantaxplor.fcm' + + INTEGER NNAT,METHOD + DOUBLE PRECISION XPX(1000),XPY(1000),XPZ(1000), + & XOBS(1000),XTOLPROT(1000),ERWPROT(1000) + DOUBLE PRECISION XPMX(1000),XPMY(1000),XPMZ(1000) + CHARACTER NOMEFILE1*132 + INTEGER ERRORE + LOGICAL LOCFLG + + PRINT*,' ' + PRINT*,'===---FANTALIN---===' +C ------------------------------------------------------------------ +C IHP is the number of PCS +C ------------------------------------------------------------------ + MDIPO=METHOD + PRINT*,'MDIPO is:',MDIPO + PRINT*,'ERRORE flag is:',ERRORE + NAT=NNAT + IHP=NNAT + NHP=NNAT + PRINT*,'IHP is:',IHP +C ------------------------------------------------------------------ +C NFE is the number of paramagnetic centers +C ------------------------------------------------------------------ + NFE=1 + NIA=0 +C ------------------------------------------------------------------ +C Reference system not active +C ------------------------------------------------------------------ + NSYSTEM=0 +C ------------------------------------------------------------------ +C Calculation is done on a single structure +C ------------------------------------------------------------------ + NSTR=1 +C ------------------------------------------------------------------ +C Errors are set to 0 +C ------------------------------------------------------------------ + PERC=0 +C ------------------------------------------------------------------ +C Grid setting +C ------------------------------------------------------------------ + NGRID=1 +C ------------------------------------------------------------------ +C Weights and multiplicity are set to 1 +C ------------------------------------------------------------------ + DO NIA=1,NAT + MLPROT(NIA)=1 + IF(ERRORE.EQ.0) THEN + WPROT(NIA)=1 + ELSE + WPROT(NIA)=ERWPROT(NIA) + END IF + END DO + + OPEN(322,FILE='FANTALIN.LOG') +C ------------------------------------------------------------------ +C Vectors are passed from Xplor-NIH core +C ------------------------------------------------------------------ + DO NIA=1,NAT + OBS(NIA)=XOBS(NIA) +C ------------------------------------------------------------------ +C The use of tolerance in fitting depends on this flag +C ------------------------------------------------------------------ + IF (LOCFLG) THEN + TOLPROT(NIA)=XTOLPROT(NIA) + ELSE + TOLPROT(NIA)=0 + END IF + CX(NSTR,NIA)=XPX(NIA) + CY(NSTR,NIA)=XPY(NIA) + CZ(NSTR,NIA)=XPZ(NIA) + WRITE (322,123) NIA,OBS(NIA),CX(NSTR,NIA),CY(NSTR,NIA), + & CZ(NSTR,NIA),TOLPROT(NIA),WPROT(NIA) + WRITE (322,55) NRESIDFA(NIA),NRESFA(NIA),NTYPEFA(NIA) +123 FORMAT (2X,I3,2X,F8.3,2X,2X,F8.3,2X,F8.3,2X,F8.3,2X,F8.3,F8.3) + 55 FORMAT (1X,'RESID',1X,A4,1X,'RES',1X,A4,1X,'TYPE',1X,A4) +C ------------------------------------------------------------------ +C Only one paramagnetic center is passed +C ------------------------------------------------------------------ + IF (METHOD.EQ.1) THEN + FX(NSTR,1)=XPMX(NIA) + FY(NSTR,1)=XPMY(NIA) + FZ(NSTR,1)=XPMZ(NIA) + ELSE + FX(NSTR,NIA)=XPMX(NIA) + FY(NSTR,NIA)=XPMY(NIA) + FZ(NSTR,NIA)=XPMZ(NIA) + END IF + WRITE (322,124) FX(NSTR,1),FY(NSTR,1),FZ(NSTR,1) +124 FORMAT (2X,F8.3,2X,F8.3,2X,F8.3) + END DO + CLOSE (322) + FILENAME3=NOMEFILE1 + PRINT*,'NOMEFILE1 is:',NOMEFILE1 + OPEN(744,FILE=FILENAME3) + + CALL SIMPLEXRUN() + CALL SCAMB(CHIPARA,CHIPERP) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE SCAMB(A,B) +C ------------------------------------------------------------------ +C The value of A2 is always negative +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + DOUBLE PRECISION A,B,APARA,APERP + COMMON /FRUN/APERP,APARA + + APARA=A + APERP=-B + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE SIMPLEXRUN() +C ------------------------------------------------------------------ +C Calculates the susceptibility tensor of paramagnetic system +C which mimimize the experimental-to-calculated squared shift +C difference +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + INCLUDE 'sup2.fcm' + + DIMENSION SIMP(MP,NP),EXTR(NP),VAL(MP) + DIMENSION ZRR(3,3) + + OLDRESID=1.E+9 + RESID=2.E+9 + IVIOLATION=2000000 + IOLDVIO=1000000 +C ------------------------------------------------------------------ +C Initial conditions grid +C ------------------------------------------------------------------ + DELTA=3.14159/DBLE(NGRID+1) + DO 1000 L=1,NGRID + DO 999 NK=1,NFE + SIMP(1,(NK-1)*MPAR+1)=1000.*(1.-2*RAND()) + SIMP(1,(NK-1)*MPAR+2)=1000.*(1.-2*RAND()) + SIMP(1,(NK-1)*MPAR+3)=1000.*(1.-2*RAND()) + SIMP(1,(NK-1)*MPAR+4)=1000.*(1.-2*RAND()) + SIMP(1,(NK-1)*MPAR+5)=1000.*(1.-2*RAND()) +999 CONTINUE + DO 10 I=1,MP + DO 10 J=1,NP + SIMP(I,J)=SIMP(1,J) + IF (I.EQ.1) THEN + SIMP(1,J)=SIMP(1,J)*0.8 + ENDIF + IF (I.EQ.J+1) THEN + SIMP(I,J)=SIMP(I,J)/0.8*1.2 + ENDIF +10 CONTINUE + DO 11 I=1,MP + DO 12 J=1,NP +12 EXTR(J)=SIMP(I,J) + VAL(I)=CASH(EXTR) +11 CONTINUE + CALL SIMPCALC(SIMP,VAL) +C ------------------------------------------------------------------ +C Criterion of minimization of optimal solution +C ------------------------------------------------------------------ + IF (OLDRESID.GE.RESID) THEN + OLDRESID=RESID + IOLDVIO=IVIOLATION + DO 13 I=1,NFE + OPTPHI(I)=PHI(I) + OPTTETA(I)=TETA(I) + OPTOMEGA(I)=OMEGA(I) + OPTA1(I)=A1DIP(I) + OPTA2(I)=A2DIP(I) +13 CONTINUE + ENDIF +1000 CONTINUE + + DO 998 NK=1,NFE + EXTR((NK-1)*MPAR+1)=OPTPHI(NK) + EXTR((NK-1)*MPAR+2)=OPTTETA(NK) + EXTR((NK-1)*MPAR+3)=OPTOMEGA(NK) + EXTR((NK-1)*MPAR+4)=OPTA1(NK) + EXTR((NK-1)*MPAR+5)=OPTA2(NK) +998 CONTINUE + + PRINT*,' ' + WRITE(*,FMT='(1X,A)') ' ***** BEST SOLUTION *****' + PRINT*,' ' + WRITE(*,20) IOLDVIO,OLDRESID + WRITE(*,21) (OPTPHI(I),OPTTETA(I),OPTOMEGA(I),I=1,NFE) + WRITE(*,22) (OPTA1(I),OPTA2(I),I=1,NFE) + + CALL DIAG(ZRR,OPTPHI(1),OPTTETA(1),OPTOMEGA(1),OPTA1(1), + & OPTA2(1),CHIPARA,CHIPERP) + + WRITE(*,*) 'A1 = ',CHIPARA,', ',CHIPARA*12*3.1416/1E4,'E-32' + WRITE(*,*) 'A2 = ',-CHIPERP,', ',-CHIPERP*12*3.1416/1E4,'E-32' + WRITE(744,*) 'A1 = ',CHIPARA,', ',CHIPARA*12*3.1416/1E4,'E-32' + WRITE(744,*) 'A2 = ',-CHIPERP,', ',-CHIPERP*12*3.1416/1E4,'E-32' + + X4=ZRR(1,1) + Y4=ZRR(1,2) + Z4=ZRR(1,3) + X4=ZRR(1,1) + Y4=ZRR(2,1) + Z4=ZRR(3,1) + X4=ZRR(2,1) + Y4=ZRR(2,2) + Z4=ZRR(2,3) + X4=ZRR(1,2) + Y4=ZRR(2,2) + Z4=ZRR(3,2) + X4=ZRR(3,1) + Y4=ZRR(3,2) + Z4=ZRR(3,3) + X4=ZRR(1,3) + Y4=ZRR(2,3) + Z4=ZRR(3,3) + DO 777 I=1,NFE + P=OPTPHI(I) + T=OPTTETA(I) + O=OPTOMEGA(I) + ASXX(I,1)=COS(P)*COS(O) + ASXX(I,2)=SIN(P)*COS(O) + ASXX(I,3)=SIN(O) + ASYY(I,1)=-COS(T)*SIN(P)-SIN(O)*COS(P)*SIN(T) + ASYY(I,2)=COS(T)*COS(P)-SIN(O)*SIN(P)*SIN(T) + ASYY(I,3)=SIN(T)*COS(O) + ASZZ(I,1)=SIN(T)*SIN(P)-SIN(O)*COS(P)*COS(T) + ASZZ(I,2)=-SIN(T)*COS(P)-SIN(O)*SIN(P)*COS(T) + ASZZ(I,3)=COS(T)*COS(O) +777 CONTINUE + + CALL PRRES + +20 FORMAT (1X,' VIOLATIONS ',I6,1X,' RESIDUALS ',F10.3) +21 FORMAT (1X,'C1= ',F15.3,1X,'C2= ',F15.3,1X,'C3= ',F15.3) +22 FORMAT (1X,'C4= ',F15.3,1X,'C5= ',F15.3) + + ERR=CASH(EXTR) + CALL DISPSHIFT + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE SIMPCALC(P,Y) +C ------------------------------------------------------------------ +C Simplex calculation +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + INCLUDE 'sup2.fcm' + INCLUDE 'fantaxplor.fcm' + + PARAMETER (AL=1.,BE=.5,GA=2.) + DIMENSION P(MP,NP),Y(MP),PR(MP),PT(MP),PE(MP), + & PC(MP),PK(MP),PBAR(MP),OLDP(NP),EXTR(MP) +C ------------------------------------------------------------------ +C Algorithm initialization +C ------------------------------------------------------------------ + ITER=0 + RNDOM=0 +1 PERM=1 + DO 1000 WHILE (PERM.GT.0) + PERM=0 + DO 1100 I=1,NP + IF (Y(I).GT.Y(I+1)) THEN + APP=Y(I+1) + Y(I+1)=Y(I) + Y(I)=APP + DO 1200 J=1,NP + APP=P(I+1,J) + P(I+1,J)=P(I,J) + P(I,J)=APP +1200 CONTINUE + PERM=PERM+1 + ENDIF +1100 CONTINUE +1000 CONTINUE + + AVERAGE=0 + ASTDEV=0 + DO 100 I=1,MP +100 AVERAGE=AVERAGE+Y(I) + AVERAGE=AVERAGE/MP + DO 101 I=1,MP +101 ASTDEV=ASTDEV+(Y(I)-AVERAGE)**2. + ERR=DSQRT(ASTDEV/MP) + RESID=Y(1) + IF (ITER.GT.39999.OR.ERR.LT.1.D-9) THEN + WRITE(*,200) ITER,IVIOLATION,RESID + WRITE (*,201) + & (P(1,(I-1)*MPAR+1), + & P(1,(I-1)*MPAR+2), + & P(1,(I-1)*MPAR+3), + & P(1,(I-1)*MPAR+4),P(1,(I-1)*MPAR+5),I=1,NFE) +200 FORMAT (1X,'ITERS ',I6,1X,' VIOLATIONS ',I6, + & 1X,'RESIDUALS= ',F10.3) +201 FORMAT (1X,'C1= ',F9.3,1X,'C2= ',F9.3,1X,'C3= ',F9.3, + & 1X,'C4= ',F9.3,1X,'C5= ',F9.3) + DO 997 NK=1,NFE + PHI(NK)=P(1,(NK-1)*MPAR+1) + TETA(NK)=P(1,(NK-1)*MPAR+2) + OMEGA(NK)=P(1,(NK-1)*MPAR+3) + A1DIP(NK)=P(1,(NK-1)*MPAR+4) + A2DIP(NK)=P(1,(NK-1)*MPAR+5) +997 CONTINUE + RETURN + ENDIF + + ITER=ITER+1 + DO 12 J=1,NP +12 PBAR(J)=0 + DO 13 I=1,NP + DO 13 J=1,NP +13 PBAR(J)=PBAR(J)+P(I,J) + DO 14 J=1,NP + PBAR(J)=PBAR(J)/NP + PR(J)=(1.+AL)*PBAR(J)-AL*P(MP,J) +14 CONTINUE + YR=CASH(PR) + DO 15 J=1,NP + PK(J)=PR(J) +15 CONTINUE + YK=YR + IF (YR.LT.Y(NP)) THEN + IF (YR.LT.Y(1)) THEN + DO 16 J=1,NP + PE(J)=GA*PR(J)+(1.-GA)*PBAR(J) +16 CONTINUE + YE=CASH(PE) + ITER=ITER+1 + IF (YE.LT.Y(1)) THEN + DO 17 J=1,NP + PK(J)=PE(J) +17 CONTINUE + YK=YE + ENDIF + ENDIF + ELSE + DO 18 J=1,NP + PT(J)=P(MP,J) +18 CONTINUE + YT=Y(MP) + IF (YR.LT.YT) THEN + DO 19 J=1,NP + PT(J)=PR(J) +19 CONTINUE + YT=YR + ENDIF + DO 20 J=1,NP + PC(J)=BE*PT(J)+(1.-BE)*PBAR(J) +20 CONTINUE + YC=CASH(PC) + ITER=ITER+1 + IF (YC.LT.Y(NP)) THEN + DO 21 J=1,NP + PK(J)=PC(J) +21 CONTINUE + YK=YC + ELSE + DO 22 I=1,NP + OLDP(I)=P(2,I) +22 CONTINUE + DO 24 I=2,NP + DO 23 J=1,NP + P(I,J)=.5*(P(1,J)+P(I,J)) +23 CONTINUE + PR(I)=P(I,I) + Y(I)=CASH(PR) +24 CONTINUE + ITER=ITER+NP-1 + DO 25 I=1,NP + OLDP(I)=ABS(OLDP(I)-P(2,I)) +25 CONTINUE + OLDPOPPA=0 + DO 251 I=1,NP-1 + OLDPOPPA=OLDPOPPA+OLDP(I) +251 CONTINUE + IF (ABS(OLDPOPPA).LT.1.D-8.AND.RNDOM.LT.10) THEN + RNDOM=RNDOM+1 + WRITE(*,FMT='(BN,A,F16.6)') 'RANDOM ',Y(1) + DO 2001 IL=1,MP + DO 2002 IM=1,NP + RMR=1+RAND() + P(IL,IM)=P(IL,IM)*RMR +2002 CONTINUE +2001 CONTINUE + DO 2011 IL=1,MP + DO 2012 JM=1,NP +2012 EXTR(JM)=P(IL,JM) + Y(IL)=CASH(EXTR) +2011 CONTINUE + GOTO 1 + ENDIF + IF (ABS(OLDPOPPA).LT.1.D-9) THEN + ITER=50000 + ENDIF + DO 26 J=1,NP + PK(J)=0.5*(P(1,J)+P(MP,J)) +26 CONTINUE + YK=CASH(PK) + ITER=ITER+1 + ENDIF + ENDIF + DO 27 J=1,NP + P(MP,J)=PK(J) +27 CONTINUE + Y(MP)=YK + GOTO 1 + END +C ------------------------------------------------------------------ + FUNCTION CASH(VETT) +C ------------------------------------------------------------------ +C Calculation of deviations and violations +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + INCLUDE 'sup2.fcm' + + DIMENSION VETT(NP) + INTEGER M1 + + IVIOLATION=0 + TMP1=0 + I=1 + DO WHILE (I.LE.NHP*NSTR) + SHIFT(I)=0.0 + I=I+1 + ENDDO + + DO 1 N=1,NSTR + DO 2 M=1,NFE + P=VETT((M-1)*MPAR+1) + T=VETT((M-1)*MPAR+2) + O=VETT((M-1)*MPAR+3) + A1D=VETT((M-1)*MPAR+4) + A2D=VETT((M-1)*MPAR+5) + TMP2=0.0 + I=1 + IHP=(N-1)*NHP + DO 10 WHILE (I.LE.NHP) + NPROTML=MLPROT(IHP+I) + IF (MDIPO.EQ.0) THEN + M1=I + ELSE + M1=M + END IF + IF (NPROTML.GT.1) THEN + COSTMULT=1/DBLE(NPROTML) + APP=0.0 + DO 11 L=1,NPROTML + IF (PRIMO.EQ.0) THEN + X4=CX(N,I)-FX(N,M1) + Y4=CY(N,I)-FY(N,M1) + Z4=CZ(N,I)-FZ(N,M1) + R(N,I)=SQRT(X4**2+Y4**2+Z4**2) + TETA1(N,I)=ACOS(Z4/R(N,I)) + PHI1(N,I)=ACOS(X4/(R(N,I)*SIN(TETA1(N,I)))) + IF (Y4.LT.0) PHI1(N,I)=2*3.1416-PHI1(N,I) + ENDIF + G1=(1.-3*COS(TETA1(N,I))**2)/R(N,I)**3 + G2=(SIN(TETA1(N,I))**2*COS(2*PHI1(N,I)))/ + & R(N,I)**3 + G3=2*(SIN(TETA1(N,I))**2*SIN(2*PHI1(N,I)))/ + & R(N,I)**3 + G4=2*(SIN(2*TETA1(N,I))*COS(PHI1(N,I)))/ + & R(N,I)**3 + G5=2*(SIN(2*TETA1(N,I))*SIN(PHI1(N,I)))/ + & R(N,I)**3 + APP=APP+3*(P*G1+T*G2+O*G3+A1D*G4+A2D*G5)/2. + I=I+1 +11 CONTINUE + DO 12 L=1,NPROTML + SHIFT(IHP+I-L)=SHIFT(IHP+I-L)+APP*COSTMULT +12 CONTINUE + ELSE + IF (PRIMO.EQ.0) THEN + X4=CX(N,I)-FX(N,M1) + Y4=CY(N,I)-FY(N,M1) + Z4=CZ(N,I)-FZ(N,M1) + R(N,I)=SQRT(X4**2+Y4**2+Z4**2) + TETA1(N,I)=ACOS(Z4/R(N,I)) + PHI1(N,I)=ACOS(X4/(R(N,I)*SIN(TETA1(N,I)))) + IF (Y4.LT.0) PHI1(N,I)=2*3.1416-PHI1(N,I) + ENDIF + G1=(1.-3*COS(TETA1(N,I))**2)/R(N,I)**3 + G2=(SIN(TETA1(N,I))**2*COS(2*PHI1(N,I)))/R(N,I)**3 + G3=2*(SIN(TETA1(N,I))**2*SIN(2*PHI1(N,I)))/R(N,I)**3 + G4=2*(SIN(2*TETA1(N,I))*COS(PHI1(N,I)))/R(N,I)**3 + G5=2*(SIN(2*TETA1(N,I))*SIN(PHI1(N,I)))/R(N,I)**3 + SHIFT(IHP+I)=SHIFT(IHP+I)+ + & 3*(P*G1+T*G2+O*G3+A1D*G4+A2D*G5)/2. + I=I+1 + ENDIF +10 CONTINUE +2 CONTINUE + DO 3 I=1,NHP + TMP2=ABS(SHIFT(IHP+I)-OBS(IHP+I))-TOLPROT(IHP+I) + IF (TMP2.GT.0.0) THEN + IVIOLATION=IVIOLATION+1 + TMP1=TMP1+TMP2**2*WPROT(IHP+I)/DBLE(MLPROT(IHP+I)) + ENDIF +3 CONTINUE +1 CONTINUE + IF (PRIMO.EQ.0) PRIMO=1 + + CASH=TMP1 + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE DISPSHIFT +C ------------------------------------------------------------------ +C Display of violations +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + INCLUDE 'sup2.fcm' + INCLUDE 'fantaxplor.fcm' + + CHARACTER FILESTR*21,SS*40 + DIMENSION VV(MAXSTR,MAXOS),VMAX(MAXOS),VMIN(MAXOS) + DIMENSION NNMAX(MAXOS),NNMIN(MAXOS) + + RM=0 + TMP1=0 + SMED=0 + STDEV=0 + DO 2 I=1,NSTR + DO 2 J=1,NHP + TMP2=ABS(OBS(J+(I-1)*NHP)-SHIFT(J+(I-1)*NHP))- + & TOLPROT(J+(I-1)*NHP) + IF (TMP2.GT.0.0) THEN + IF (SHIFT(J+(I-1)*NHP).GT.OBS(J+(I-1)*NHP)) THEN + VV(I,J)=TMP2 + ELSE + VV(I,J)=-TMP2 + ENDIF + ELSE + VV(I,J)=0. + ENDIF +2 CONTINUE + DO 44 I=1,NHP + VMAX(I)=0. + VMIN(I)=1000. +44 CONTINUE + DO 4 J=1,NHP + DO 4 I=1,NSTR + IF (ABS(VV(I,J)).GT.ABS(VMAX(J))) THEN + VMAX(J)=VV(I,J) + ENDIF + IF ((ABS(VV(I,J)).LT.ABS(VMIN(J))).AND.(VV(I,J).NE.0.)) THEN + VMIN(J)=VV(I,J) + ENDIF +4 CONTINUE + DO 7 J=1,NHP + DO 8 I=1,NSTR + IF (ABS(VV(I,J)).GT.0.) THEN + SS(I:I+1)="*" + ELSE + SS(I:I+1)=" " + ENDIF +8 CONTINUE + IF (VMIN(J).EQ.1000.) THEN + VMIN(J)=0. + ENDIF + +7 CONTINUE + + WRITE(744,'(1X)') + WRITE(744,'(1X)') + WRITE(744,'(1X)') + WRITE(744,'(1X,A)') + & ' PROTONS OBS. CALC. ERR.^2' + + FILESTR='USELESS.FILE' + NIA=0 + REWIND(1) + DO 3 J=1,NHP*NSTR + TMP2=ABS(OBS(J)-SHIFT(J))-TOLPROT(J) + NIA=NIA+1 + IF (TMP2.GT.0.0) THEN + IVIOLATION=IVIOLATION+1 + RM=RM+1 + SMED=SMED+TMP2**2*WPROT(J)/MLPROT(J) + TMP1=TMP2**2*WPROT(J)/MLPROT(J) + WRITE (744,'(1X,A4,1X,A4,1X,A4,1X,F7.3,1X,F7.3,1X,F10.5)') + & NRESIDFA(NIA),NRESFA(NIA),NTYPEFA(NIA), + & OBS(J),SHIFT(J),TMP1 +15 FORMAT(1X,A,1X,F7.3,1X,F7.3,1X,F7.3) + ELSE + WRITE (744,'(1X,A4,1X,A4,1X,A4,1X,F7.3,1X,F7.3)') + & NRESIDFA(NIA),NRESFA(NIA),NTYPEFA(NIA), + & OBS(J),SHIFT(J) + ENDIF +3 CONTINUE + SMED=SMED/RM + + DO 41 J=1,NHP*NSTR + TMP2=ABS(OBS(J)-SHIFT(J))-TOLPROT(J) + IF (TMP2.GT.0.0) THEN + STDEV=STDEV+(TMP2**2-MED)**2*WPROT(J)/MLPROT(J) + ENDIF +41 CONTINUE + STDEV=STDEV/RM + + PRINT*,' ' + WRITE(744,'(A)') ' ' + WRITE(*,11) 'Mean value of error =',SMED +11 FORMAT(1X,A,1X,F10.6) + WRITE(744,12) 'Mean value of error =',SMED +12 FORMAT(1X,A,1X,F10.6) + WRITE(*,13) 'Standard deviation =',STDEV +13 FORMAT(1X,A,1X,F10.6) + WRITE(744,14) 'Standard deviation =',STDEV +14 FORMAT(1X,A,1X,F10.6) + CLOSE(744) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE PRRES +C ------------------------------------------------------------------ +C SAXX, SAYY, SAZZ : Direction cosines of rotation matrices +C SAXSYS, SAYSYS, SAZSYS : Direction cosines of reference systems +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + INCLUDE 'sup2.fcm' + INCLUDE 'fantaxplor.fcm' + + DIMENSION SAXX(3),SAYY(3),SAZZ(3), + & APPX(3),APPY(3),APPZ(3), + & SAXSYS(3),SAYSYS(3),SAZSYS(3) + CHARACTER NAMATOM*4,NAMRESS*3 + PI=3.1415927 + +99 I=0 + DO 10 N=1,NFE + A1=OPTA1(N) + A2=OPTA2(N) + IF (NSYSTEM.GT.1) THEN + DO 1 I=1,3 + SAXSYS(I)=AXSYS(N,I) + SAYSYS(I)=AYSYS(N,I) + SAZSYS(I)=AZSYS(N,I) +1 CONTINUE + ELSE + DO 2 I=1,3 + SAXSYS(I)=AXSYS(1,I) + SAYSYS(I)=AYSYS(1,I) + SAZSYS(I)=AZSYS(1,I) +2 CONTINUE + ENDIF + DO 3 I=1,3 + SAXX(I)=ASXX(N,I) + SAYY(I)=ASYY(N,I) + SAZZ(I)=ASZZ(N,I) +3 CONTINUE + SDZX=SAXX(2)*SAYY(3)-SAXX(3)*SAYY(2) + IF (SDZX/SAZZ(1).GT.0.0) THEN +C ------------------------------------------------------------------ +C Right-handed +C ------------------------------------------------------------------ + PRINT*,' ' + ELSE +C ------------------------------------------------------------------ +C Left-handed +C ------------------------------------------------------------------ + PRINT *,' ' + DO 4 I=1,3 + APPX(I)=SAXX(I) + APPY(I)=SAYY(I) +4 CONTINUE + DO 5 I=1,3 + SAXX(I)=APPY(I) + SAYY(I)=APPX(I) +5 CONTINUE + A2=-A2 + ENDIF + PRINT*, + & '***********************************************************' + WRITE(744,'(A)') ' ' + WRITE(744,'(A)') + & '**********************************************************' + NUMBERATOM=NUMBERATOM+1 + NAMATOM=' CEN' + NAMRESS='PMC' + NORES=NORES+N +10 CONTINUE +22 FORMAT (1X,'A1= ',F9.3,2X,A1,F7.3,A,2X,'A2= ',F9.3,2X,A1,F7.3,A) +23 FORMAT (1X,'ARCCOS X^XT = ',F8.3,4X,'ARCCOS Y^YT = ',F8.3, + & 4X,'ARCCOS Z^ZT = ',F8.3) + + END +C ------------------------------------------------------------------ + SUBROUTINE PANG(N,NR) +C ------------------------------------------------------------------ +C Angle permutations +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + INCLUDE 'sup2.fcm' + INCLUDE 'fantaxplor.fcm' + PI=3.1415927 + + IF (N.EQ.0) THEN +C ------------------------------------------------------------------ +C No permutation +C ------------------------------------------------------------------ + IF (ASXX(NR,3).EQ.0.0) THEN + O=0.0 + P=ASIN(ASXX(NR,2)) + T=ASIN(DSQRT(ASZZ(NR,1)**2+ASZZ(NR,2)**2)) + ELSE IF (ASXX(NR,3).EQ.1.0) THEN + O=PI/2 + P=ASIN(DSQRT(ASYY(NR,1)**2+ASZZ(NR,1)**2)) + T=ASIN(DSQRT(ASZZ(NR,1)**2+ASZZ(NR,2)**2)) + ELSE IF (ASXX(NR,3).EQ.-1.0) THEN + O=-PI/2 + P=ASIN(DSQRT(ASYY(NR,1)**2+ASZZ(NR,1)**2)) + T=ASIN(DSQRT(ASZZ(NR,1)**2+ASZZ(NR,2)**2)) + ELSE + O=ASIN(ASXX(NR,3)) + P=ASIN(DSQRT((ASYY(NR,1)**2+ASZZ(NR,1)**2-ASXX(NR,3)**2)/ + & (1.-ASXX(NR,3)**2))) + T=ASIN(DSQRT((ASZZ(NR,1)**2+ASZZ(NR,2)**2-ASXX(NR,3)**2)/ + & (1.-ASXX(NR,3)**2))) + ENDIF + ELSE IF (N.EQ.1) THEN +C ------------------------------------------------------------------ +C Permutation 1 +C ------------------------------------------------------------------ + IF (ASZZ(NR,3).EQ.0.0) THEN + O=0.0 + P=ASIN(ASZZ(NR,2)) + T=ASIN(DSQRT(ASYY(NR,1)**2+ASYY(NR,2)**2)) + ELSE IF (ASZZ(NR,3).EQ.1.0) THEN + O=PI/2 + P=ASIN(DSQRT(ASXX(NR,1)**2+ASYY(NR,1)**2)) + T=ASIN(DSQRT(ASYY(NR,1)**2+ASYY(NR,2)**2)) + ELSE IF (ASZZ(NR,3).EQ.-1.0) THEN + O=-PI/2 + P=ASIN(DSQRT(ASXX(NR,1)**2+ASYY(NR,1)**2)) + T=ASIN(DSQRT(ASYY(NR,1)**2+ASYY(NR,2)**2)) + ELSE + O=ASIN(ASZZ(NR,3)) + P=DASIN(DSQRT((ASXX(NR,1)**2+ASYY(NR,1)**2-ASZZ(NR,3)**2)/ + & (1.-ASZZ(NR,3)**2))) + T=ASIN(DSQRT((ASYY(NR,1)**2+ASYY(NR,2)**2-ASZZ(NR,3)**2)/ + & (1.-ASZZ(NR,3)**2))) + ENDIF + ELSE +C ------------------------------------------------------------------ +C Permutation 2 +C ------------------------------------------------------------------ + IF (ASYY(NR,3).EQ.0.0) THEN + O=0.0 + P=ASIN(ASYY(NR,2)) + T=ASIN(DSQRT(ASXX(NR,1)**2+ASXX(NR,2)**2)) + ELSE IF (ASYY(NR,3).EQ.1.0) THEN + O=PI/2 + P=ASIN(DSQRT(ASZZ(NR,1)**2+ASXX(NR,1)**2)) + T=ASIN(DSQRT(ASXX(NR,1)**2+ASXX(NR,2)**2)) + ELSE IF (ASYY(NR,3).EQ.-1.0) THEN + O=-PI/2 + P=ASIN(DSQRT(ASZZ(NR,1)**2+ASXX(NR,1)**2)) + T=ASIN(DSQRT(ASXX(NR,1)**2+ASXX(NR,2)**2)) + ELSE + O=ASIN(ASYY(NR,3)) + P=ASIN(DSQRT((ASZZ(NR,1)**2+ASXX(NR,1)**2-ASYY(NR,3)**2)/ + & (1.-ASYY(NR,3)**2))) + T=ASIN(DSQRT((ASXX(NR,1)**2+ASXX(NR,2)**2-ASYY(NR,3)**2)/ + & (1.-ASYY(NR,3)**2))) + ENDIF + ENDIF + + SXX=COS(P)*COS(O) + SXY=SIN(P)*COS(O) + SXZ=SIN(O) + SYX=-COS(T)*SIN(P)-SIN(O)*COS(P)*SIN(T) + SYY=COS(T)*COS(P)-SIN(O)*SIN(P)*SIN(T) + SYZ=SIN(T)*COS(O) + SZX=SIN(T)*SIN(P)-SIN(O)*COS(P)*COS(T) + SZY=-SIN(T)*COS(P)-SIN(O)*SIN(P)*COS(T) + SZZ=COS(T)*COS(O) + VZX=SXY*SYZ-SXZ*SYY + VZY=SXZ*SYX-SXX*SYZ + VZZ=SXX*SYY-SXY*SYX + CDIR=SZX*VZX+SZY*VZY+SZZ*VZZ + IF (CDIR.LT.0.0) THEN + P=-P + PRINT*,'Changing PHI' + ENDIF + SXX=COS(P)*COS(O) + SXY=SIN(P)*COS(O) + SXZ=SIN(O) + SYX=-COS(T)*SIN(P)-SIN(O)*COS(P)*SIN(T) + SYY=COS(T)*COS(P)-SIN(O)*SIN(P)*SIN(T) + SYZ=SIN(T)*COS(O) + SZX=SIN(T)*SIN(P)-SIN(O)*COS(P)*COS(T) + SZY=-SIN(T)*COS(P)-SIN(O)*SIN(P)*COS(T) + SZZ=COS(T)*COS(O) + VZX=SXY*SYZ-SXZ*SYY + VZY=SXZ*SYX-SXX*SYZ + VZZ=SXX*SYY-SXY*SYX + CDIR=SZX*VZX+SZY*VZY+SZZ*VZZ + IF (CDIR.LT.0.0) THEN + T=-T + PRINT*,'Changing THETA' + ENDIF + WRITE (*,21) P,T,O +21 FORMAT(1X,'PHI= ',F7.3,4X,'THETA= ',F7.3,4X,'OMEGA= ',F7.3) + + END +C ------------------------------------------------------------------ + SUBROUTINE DIAG(ZRR,R1,R2,R3,R4,R5,CHIPARA,CHIPERP) +C ------------------------------------------------------------------ +C Handles matrix diagonalization +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT REAL*8(A-H,O-Z) + DIMENSION WR(3) + DIMENSION ARR(3,3),ARI(3,3) + DIMENSION ZRR(3,3),ZRI(3,3) + DIMENSION WK1(3),WK2(3),WK3(3) +C ------------------------------------------------------------------ +C Energy matrix for hyperfine coupling +C ------------------------------------------------------------------ + ARR(1,1)=0 + ARR(2,2)=-R2 + ARR(3,3)=(ARR(2,2)-3*R1)/2 + ARR(1,2)=R3 + ARR(1,3)=R4 + ARR(2,3)=R5 + ARR(2,1)=ARR(1,2) + ARR(3,2)=ARR(2,3) + ARR(3,1)=ARR(1,3) + NMX=3 + MBRANC=3 +C ------------------------------------------------------------------ +C Energy matrix diagonalization +C ------------------------------------------------------------------ + DO 45 I=1,NMX + DO 45 J=1,NMX + ARI(I,J)=0 +45 CONTINUE + + CALL F02AXF(ARR,MBRANC,ARI,MBRANC,NMX,WR,ZRR,MBRANC,ZRI,MBRANC, + & WK1,WK2,WK3,0) + + CHIPARA=WR(3)-(WR(2)+WR(1))/2 + CHIPERP=WR(2)-WR(1) + IF (DABS(CHIPERP).GT.2./3.*DABS(CHIPARA)) THEN + CHIPARA=DABS(WR(1)-(WR(3)+WR(2))/2) + CHIPERP=DABS(WR(3)-WR(2)) + END IF + WRITE(6,*)'Eigenvalues:',WR(1),WR(2),WR(3) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE F01BCF(N,TOL,Z,IZ,W,IW,D,E,C,S) +C ------------------------------------------------------------------ +C MARK 3 RELEASE. NAG COPYRIGHT 1972. +C MARK 4 REVISED. +C MARK 4.5 REVISED. +C MARK 5C REVISED. +C MARK 6B REVISED IER-125 IER-127 (Jul 1978). +C MARK 11 REVISED VECTORISATION (Jan 1984). +C MARK 11.5 (F77) REVISED (Sept 1985). +C +C Reduces a complex Hermitian matrix to real tridiagonal form from +C which the eigenvalues and eigenvectors can be found using +C subroutine F02AYF. The Hermitian matrix A=A(1) is reduced to the +C tridiagonal matrix A(N-1) by (N-2) unitary transformations. +C The Householder reduction itself does not give a real tridiagonal +C matrix, because the off-diagonal elements are complex. They are +C subsequently made real by a diagonal transformation. +C +C April 1st. 1972 +C ------------------------------------------------------------------ +C Scalar arguments +C ------------------------------------------------------------------ + DOUBLE PRECISION TOL + INTEGER IW, IZ, N +C ------------------------------------------------------------------ +C Array arguments +C ------------------------------------------------------------------ + DOUBLE PRECISION C(N), D(N), E(N), S(N), W(IW,N), Z(IZ,N) +C ------------------------------------------------------------------ +C Local scalars +C ------------------------------------------------------------------ + DOUBLE PRECISION CO, F, FI, FR, G, GI, GR, H, HH, R, SI + INTEGER I, II, J, K, L +C ------------------------------------------------------------------ +C External subroutines +C ------------------------------------------------------------------ + EXTERNAL F01BCY, F01BCZ +C ------------------------------------------------------------------ +C Intrinsic functions +C ------------------------------------------------------------------ + INTRINSIC ABS, SQRT +C ------------------------------------------------------------------ +C Executable statements +C ------------------------------------------------------------------ + DO 20 I = 1, N + D(I) = Z(N,I) + E(I) = -W(N,I) +20 CONTINUE + IF (N.EQ.1) GO TO 540 + DO 360 II = 2, N + I = N - II + 2 + L = I - 2 + G = 0.0D0 + FR = D(I-1) + FI = E(I-1) + IF (L.EQ.0) GO TO 60 + DO 40 K = 1, L + G = G + D(K)*D(K) + E(K)*E(K) +40 CONTINUE +60 H = G + FR*FR + FI*FI +C ------------------------------------------------------------------ +C L is now I-1 +C ------------------------------------------------------------------ + L = L + 1 + IF (ABS(FR)+ABS(FI).NE.0.0D0) GO TO 80 + R = 0.0D0 + CO = 1.0D0 + C(I) = 1.0D0 + SI = 0.0D0 + S(I) = 0.0D0 + GO TO 140 +80 IF (ABS(FR).LT.ABS(FI)) GO TO 100 + R = ABS(FR)*SQRT(1.0D0+(FI/FR)**2) + GO TO 120 +100 R = ABS(FI)*SQRT(1.0D0+(FR/FI)**2) +120 SI = FI/R + S(I) = -SI + CO = FR/R + C(I) = CO +140 IF (G.LE.TOL) GO TO 280 + G = -SQRT(H) + E(I) = G +C ------------------------------------------------------------------ +C E(I) has its final real value +C ------------------------------------------------------------------ + H = H - R*G +C ------------------------------------------------------------------ +C S*S + SR +C ------------------------------------------------------------------ + D(I-1) = (R-G)*CO + E(I-1) = (R-G)*SI + DO 160 J = 1, L + Z(J,I) = D(J) + W(J,I) = E(J) +160 CONTINUE + + CALL F01BCZ(Z,IZ,W,IW,L,D,E,C,S) +C ------------------------------------------------------------------ +C Form P +C ------------------------------------------------------------------ + DO 180 J = 1, L + C(J) = C(J)/H + S(J) = S(J)/H +180 CONTINUE + FR = 0.0D0 + DO 200 J = 1, L + FR = FR + C(J)*D(J) + S(J)*E(J) +200 CONTINUE +C ------------------------------------------------------------------ +C Form K +C ------------------------------------------------------------------ + HH = FR/(H+H) +C ------------------------------------------------------------------ +C Form Q +C ------------------------------------------------------------------ + DO 220 J = 1, L + C(J) = C(J) - HH*D(J) + S(J) = S(J) - HH*E(J) +220 CONTINUE +C ------------------------------------------------------------------ +C Now form reduced A +C ------------------------------------------------------------------ + DO 260 J = 1, L + FR = D(J) + FI = E(J) + GR = C(J) + GI = S(J) + DO 240 K = J, L + Z(K,J) = (((Z(K,J)-GR*D(K))-GI*E(K))-FR*C(K)) - FI*S(K) + W(K,J) = (((W(K,J)-GR*E(K))+GI*D(K))-FR*S(K)) + FI*C(K) +240 CONTINUE + D(J) = Z(L,J) + Z(I,J) = 0.0D0 + E(J) = -W(L,J) + W(I,J) = 0.0D0 + W(J,J) = 0.0D0 +260 CONTINUE + GO TO 340 +280 E(I) = R + H = 0.0D0 + DO 300 J = 1, L + Z(J,I) = D(J) + W(J,I) = E(J) +300 CONTINUE + DO 320 J = 1, L + Z(I,J) = 0.0D0 + D(J) = Z(I-1,J) + W(I,J) = 0.0D0 + E(J) = -W(I-1,J) +320 CONTINUE +340 D(I) = H +360 CONTINUE +C ------------------------------------------------------------------ +C We now form the product of the Householder matrices, +C overwriting on Z and W +C ------------------------------------------------------------------ + DO 500 I = 2, N + L = I - 1 + Z(N,L) = Z(L,L) + Z(L,L) = 1.0D0 + W(N,L) = E(L) + W(L,L) = 0.0D0 + H = D(I) + IF (H.EQ.0.0D0) GO TO 460 + DO 380 K = 1, L + D(K) = 0.0D0 + E(K) = 0.0D0 +380 CONTINUE + + CALL F01BCY(Z,IZ,W,IW,L,L,Z(1,I),W(1,I),D,E) + + DO 400 K = 1, L + D(K) = D(K)/H + E(K) = -E(K)/H +400 CONTINUE + DO 440 J = 1, L + DO 420 K = 1, L + Z(K,J) = Z(K,J) - Z(K,I)*D(J) + W(K,I)*E(J) + W(K,J) = W(K,J) - Z(K,I)*E(J) - W(K,I)*D(J) +420 CONTINUE +440 CONTINUE +460 DO 480 J = 1, L + Z(J,I) = 0.0D0 + W(J,I) = 0.0D0 +480 CONTINUE +500 CONTINUE + W(N,N) = E(N) + DO 520 I = 1, N + D(I) = Z(N,I) + Z(N,I) = 0.0D0 + E(I) = W(N,I) + W(N,I) = 0.0D0 +520 CONTINUE +540 Z(N,N) = 1.0D0 + W(N,N) = 0.0D0 + E(1) = 0.0D0 +C ------------------------------------------------------------------ +C Now we multiply by the COSTHETA + I SINTHETA column factors +C ------------------------------------------------------------------ + CO = 1.0D0 + SI = 0.0D0 + IF (N.EQ.1) RETURN + DO 580 I = 2, N + F = CO*C(I) - SI*S(I) + SI = CO*S(I) + SI*C(I) + CO = F + DO 560 J = 1, N + F = Z(J,I)*CO - W(J,I)*SI + W(J,I) = Z(J,I)*SI + W(J,I)*CO + Z(J,I) = F +560 CONTINUE +580 CONTINUE + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE F01BCY(AR,IAR,AI,IAI,M,N,BR,BI,CR,CI) +C ------------------------------------------------------------------ +C MARK 11 RELEASE. NAG COPYRIGHT 1983. +C MARK 11.5 (F77) REVISED (SEPT 1985). +C +C Computes C=C+(A**H)*B (COMPLEX) where: +C A is rectangular M by N +C C must be distinct from B +C ------------------------------------------------------------------ +C Scalar arguments +C ------------------------------------------------------------------ + INTEGER IAI, IAR, M, N +C ------------------------------------------------------------------ +C Array arguments +C ------------------------------------------------------------------ + DOUBLE PRECISION AI(IAI,N), AR(IAR,N), BI(M), BR(M), CI(N), CR(N) +C ------------------------------------------------------------------ +C Local scalars +C ------------------------------------------------------------------ + DOUBLE PRECISION XI, XR + INTEGER I, J +C ------------------------------------------------------------------ +C Executable statements +C ------------------------------------------------------------------ + DO 40 I = 1, N + XR = CR(I) + XI = CI(I) + DO 20 J = 1, M + XR = XR + AR(J,I)*BR(J) + AI(J,I)*BI(J) + XI = XI + AR(J,I)*BI(J) - AI(J,I)*BR(J) +20 CONTINUE + CR(I) = XR + CI(I) = XI +40 CONTINUE + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE F01BCZ(AR,IAR,AI,IAI,N,BR,BI,CR,CI) +C ------------------------------------------------------------------ +C MARK 11 RELEASE. NAG COPYRIGHT 1983. +C MARK 11.5 (F77) REVISED (SEPT 1985). +C +C Computes C=A*B (COMPLEX) where: +C A is a Hermitian N-by-N matrix whose lower triangle is stored in A +C C must be distinct from B +C ------------------------------------------------------------------ +C Scalar arguments +C ------------------------------------------------------------------ + INTEGER IAI, IAR, N +C ------------------------------------------------------------------ +C Array arguments +C ------------------------------------------------------------------ + DOUBLE PRECISION AI(IAI,N), AR(IAR,N), BI(N), BR(N), CI(N), CR(N) +C ------------------------------------------------------------------ +C Local scalars +C ------------------------------------------------------------------ + DOUBLE PRECISION YI, YR + INTEGER I, IP1, J, NM1 +C ------------------------------------------------------------------ +C Executable statements +C ------------------------------------------------------------------ + DO 20 I = 1, N + CR(I) = 0.0D0 + CI(I) = 0.0D0 +20 CONTINUE + IF (N.EQ.1) GO TO 100 + NM1 = N - 1 + DO 80 I = 1, NM1 + DO 40 J = I, N + CR(J) = CR(J) + AR(J,I)*BR(I) - AI(J,I)*BI(I) + CI(J) = CI(J) + AR(J,I)*BI(I) + AI(J,I)*BR(I) +40 CONTINUE + YR = CR(I) + YI = CI(I) + IP1 = I + 1 + DO 60 J = IP1, N + YR = YR + AR(J,I)*BR(J) + AI(J,I)*BI(J) + YI = YI + AR(J,I)*BI(J) - AI(J,I)*BR(J) +60 CONTINUE + CR(I) = YR + CI(I) = YI +80 CONTINUE +100 CR(N) = CR(N) + AR(N,N)*BR(N) - AI(N,N)*BI(N) + CI(N) = CI(N) + AR(N,N)*BI(N) + AI(N,N)*BR(N) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE F02AXF(AR,IAR,AI,IAI,N,WR,VR,IVR,VI,IVI,WK1,WK2,WK3, + & IFAIL) +C ------------------------------------------------------------------ +C MARK 3 RELEASE. NAG COPYRIGHT 1972. +C MARK 4.5 REVISED. +C MARK 9 REVISED IER-327 (Sep 1981). +C MARK 11.5 (F77) REVISED (Sep 1985). +C MARK 13 REVISED. USE OF MARK 12 X02 FUNCTIONS (Apr 1988). +C +C Eigenvalues and eigenvectors of a complex Hermitian matrix +C 1st April 1972 +C ------------------------------------------------------------------ +C Parameters +C ------------------------------------------------------------------ + CHARACTER*6 SRNAME + PARAMETER (SRNAME='F02AXF') +C ------------------------------------------------------------------ +C Scalar arguments +C ------------------------------------------------------------------ + INTEGER IAI, IAR, IFAIL, IVI, IVR, N +C ------------------------------------------------------------------ +C Array arguments +C ------------------------------------------------------------------ + DOUBLE PRECISION AI(IAI,N), AR(IAR,N), VI(IVI,N), VR(IVR,N), + & WK1(N), WK2(N), WK3(N), WR(N) +C ------------------------------------------------------------------ +C Local scalars +C ------------------------------------------------------------------ + DOUBLE PRECISION A, B, MAX, SQ, SUM, XXXX + INTEGER I, ISAVE, J +C ------------------------------------------------------------------ +C Locla arrays +C ------------------------------------------------------------------ + CHARACTER*1 P01REC(1) +C ------------------------------------------------------------------ +C External functions +C ------------------------------------------------------------------ + DOUBLE PRECISION X02AJF, X02AKF + INTEGER P01ABF + EXTERNAL X02AJF, X02AKF, P01ABF +C ------------------------------------------------------------------ +C External subroutines +C ------------------------------------------------------------------ + EXTERNAL F01BCF, F02AYF +C ------------------------------------------------------------------ +C Intrinsic functions +C ------------------------------------------------------------------ + INTRINSIC SQRT +C ------------------------------------------------------------------ +C Executable statements +C ------------------------------------------------------------------ + ISAVE = IFAIL + DO 40 I = 1, N + IF (AI(I,I).NE.0.0D0) GO TO 140 + DO 20 J = 1, I + VR(I,J) = AR(I,J) + VI(I,J) = AI(I,J) +20 CONTINUE +40 CONTINUE + + CALL F01BCF(N,X02AKF()/X02AJF(),VR,IVR,VI,IVI,WR,WK1,WK2,WK3) + CALL F02AYF(N,X02AJF(),WR,WK1,VR,IVR,VI,IVI,IFAIL) + + IF (IFAIL.EQ.0) GO TO 60 + IFAIL = P01ABF(ISAVE,1,SRNAME,0,P01REC) + RETURN + +60 DO 120 I = 1, N + SUM = 0.0D0 + MAX = 0.0D0 + DO 80 J = 1, N + SQ = VR(J,I)*VR(J,I) + VI(J,I)*VI(J,I) + SUM = SUM + SQ + IF (SQ.LE.MAX) GO TO 80 + MAX = SQ + A = VR(J,I) + B = VI(J,I) +80 CONTINUE + IF (SUM.EQ.0.0D0) GO TO 120 + SUM = 1.0D0/SQRT(SUM*MAX) + DO 100 J = 1, N + SQ = SUM*(VR(J,I)*A+VI(J,I)*B) + VI(J,I) = SUM*(VI(J,I)*A-VR(J,I)*B) + VR(J,I) = SQ +100 CONTINUE +120 CONTINUE + RETURN +140 IFAIL = P01ABF(ISAVE,2,SRNAME,0,P01REC) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE F02AYF(N,EPS,D,E,Z,IZ,W,IW,IFAIL) +C ------------------------------------------------------------------ +C MARK 3 RELEASE. NAG COPYRIGHT 1972. +C MARK 4 REVISED. +C MARK 4.5 REVISED. +C MARK 9 REVISED IER-326 (SEP 1981). +C MARK 11.5 (F77) REVISED (SEP 1985). +C +C This subroutine finds the eigenvalues and eigenvectors of a +C Hermitian matrix, which has been reduced to a real tridiagonal +C matrix, T, given with its diagonal elements in the array D(N) and +C its sub-diagonal elements in the last N-1 stores of the array +C E(N), using QL transformations. The eigenvalues are overwritten +C on the diagonal elements in the array D in ascending order. The +C real and imaginary parts of the eigenvectors are formed in the +C arrays Z,W(N,N) respectively, overwriting the accumulated +C transformations as supplied by the subroutine F01BCF. The +C subroutine will fail if all eigenvalues take more than 30*N +C iterations. +C 1st April 1972 +C ------------------------------------------------------------------ +C Parameters +C ------------------------------------------------------------------ + CHARACTER*6 SRNAME + PARAMETER (SRNAME='F02AYF') +C ------------------------------------------------------------------ +C Scalar arguments +C ------------------------------------------------------------------ + DOUBLE PRECISION EPS + INTEGER IFAIL, IW, IZ, N +C ------------------------------------------------------------------ +C Array arguments +C ------------------------------------------------------------------ + DOUBLE PRECISION D(N), E(N), W(IW,N), Z(IZ,N) +C ------------------------------------------------------------------ +C Local scalars +C ------------------------------------------------------------------ + DOUBLE PRECISION B, C, F, G, H, P, R, S + INTEGER I, I1, II, ISAVE, J, K, L, M, M1 +C ------------------------------------------------------------------ +C Local arrays +C ----------------------------------------------------------------- + CHARACTER*1 P01REC(1) +C ------------------------------------------------------------------ +C External functions +C ------------------------------------------------------------------ + INTEGER P01ABF + EXTERNAL P01ABF +C ------------------------------------------------------------------ +C Intrinsic functions +C ------------------------------------------------------------------ + INTRINSIC ABS, SQRT +C ------------------------------------------------------------------ +C Executable statements +C ------------------------------------------------------------------ + ISAVE = IFAIL + IF (N.EQ.1) GO TO 40 + DO 20 I = 2, N + E(I-1) = E(I) +20 CONTINUE +40 E(N) = 0.0D0 + B = 0.0D0 + F = 0.0D0 + J = 30*N + DO 300 L = 1, N + H = EPS*(ABS(D(L))+ABS(E(L))) + IF (B.LT.H) B = H +C ------------------------------------------------------------------ +C Look for small sub-diagonal elements +C ------------------------------------------------------------------ + DO 60 M = L, N + IF (ABS(E(M)).LE.B) GO TO 80 +60 CONTINUE +80 IF (M.EQ.L) GO TO 280 +100 IF (J.LE.0) GO TO 400 + J = J - 1 +C ------------------------------------------------------------------ +C Form shift +C ------------------------------------------------------------------ + G = D(L) + H = D(L+1) - G + IF (ABS(H).GE.ABS(E(L))) GO TO 120 + P = H*0.5D0/E(L) + R = SQRT(P*P+1.0D0) + H = P + R + IF (P.LT.0.0D0) H = P - R + D(L) = E(L)/H + GO TO 140 +120 P = 2.0D0*E(L)/H + R = SQRT(P*P+1.0D0) + D(L) = E(L)*P/(1.0D0+R) +140 H = G - D(L) + I1 = L + 1 + IF (I1.GT.N) GO TO 180 + DO 160 I = I1, N + D(I) = D(I) - H +160 CONTINUE +180 F = F + H +C ------------------------------------------------------------------ +C QL transformation +C ------------------------------------------------------------------ + P = D(M) + C = 1.0D0 + S = 0.0D0 + M1 = M - 1 + DO 260 II = L, M1 + I = M1 - II + L + G = C*E(I) + H = C*P + IF (ABS(P).LT.ABS(E(I))) GO TO 200 + C = E(I)/P + R = SQRT(C*C+1.0D0) + E(I+1) = S*P*R + S = C/R + C = 1.0D0/R + GO TO 220 +200 C = P/E(I) + R = SQRT(C*C+1.0D0) + E(I+1) = S*E(I)*R + S = 1.0D0/R + C = C/R +220 P = C*D(I) - S*G + D(I+1) = H + S*(C*G+S*D(I)) +C ------------------------------------------------------------------ +C Form vector +C ------------------------------------------------------------------ + DO 240 K = 1, N + H = Z(K,I+1) + Z(K,I+1) = S*Z(K,I) + C*H + Z(K,I) = C*Z(K,I) - S*H + H = W(K,I+1) + W(K,I+1) = S*W(K,I) + C*H + W(K,I) = C*W(K,I) - S*H +240 CONTINUE +260 CONTINUE + E(L) = S*P + D(L) = C*P + IF (ABS(E(L)).GT.B) GO TO 100 +280 D(L) = D(L) + F +300 CONTINUE +C ------------------------------------------------------------------ +C Sort eigenvalues and eigenvectors +C ------------------------------------------------------------------ + DO 380 I = 1, N + K = I + P = D(I) + I1 = I + 1 + IF (I1.GT.N) GO TO 340 + DO 320 J = I1, N + IF (D(J).GE.P) GO TO 320 + K = J + P = D(J) +320 CONTINUE +340 IF (K.EQ.I) GO TO 380 + D(K) = D(I) + D(I) = P + DO 360 J = 1, N + P = Z(J,I) + Z(J,I) = Z(J,K) + Z(J,K) = P + P = W(J,I) + W(J,I) = W(J,K) + W(J,K) = P +360 CONTINUE +380 CONTINUE + RETURN +400 IFAIL = P01ABF(ISAVE,1,SRNAME,0,P01REC) + + RETURN + END +C ------------------------------------------------------------------ + INTEGER FUNCTION P01ABF(IFAIL,IERROR,SRNAME,NREC,REC) +C ------------------------------------------------------------------ +C MARK 11.5 (F77) RELEASE. NAG COPYRIGHT 1986. +C MARK 13 REVISED IER-621 (APR 1988). +C MARK 13B REVISED IER-668 (AUG 1988). +C +C P01ABF is the error-handling routine for the NAG library. +C +C P01ABF either returns the value of IERROR through the routine +C name (soft failure), or terminates execution of the program +C (hard failure). Diagnostic messages may be output. +C +C If IERROR = 0 (successful exit from the calling routine), +C the value 0 is returned through the routine name, and no +C message is ouptut. +C +C If IERROR is non-zero (abnormal exit from the calling routine), +C the action taken depends on the value of IFAIL. +C +C IFAIL = 1: Soft failure, silent exit (i.e. no messages are +C output) +C IFAIL = -1: Soft failure, noisy exit (i.e. messages are output) +C IFAIL =-13: Soft failure, noisy exit but standard messages from +C P01ABF are suppressed +C IFAIL = 0: Hard failure, noisy exit +C +C For compatibility with certain routines included before MARK 12 +C P01ABF also allows an alternative specification of IFAIL in which +C it is regarded as a decimal integer with least significant digits +C CBA. Then: +C +C A = 0: Hard failure A = 1: Soft failure +C B = 0: Silent exit B = 1: Noisy exit +C +C Except that hrad failure now always implies a noisy exit. +C +C S.Hammarling, M.P.Hooper and J.J.Du Croz, NAG Central Office. +C ------------------------------------------------------------------ +C Scalar arguments +C ------------------------------------------------------------------ + INTEGER IERROR, IFAIL, NREC + CHARACTER*(*) SRNAME +C ------------------------------------------------------------------ +C Array arguments +C ------------------------------------------------------------------ + CHARACTER*(*) REC(*) +C ------------------------------------------------------------------ +C Local scalars +C ------------------------------------------------------------------ + INTEGER I, NERR + CHARACTER*72 MESS +C ------------------------------------------------------------------ +C External subroutines +C ------------------------------------------------------------------ + EXTERNAL P01ABZ, X04AAF, X04BAF +C ------------------------------------------------------------------ +C Intrinsic functions +C ------------------------------------------------------------------ + INTRINSIC ABS, MOD +C ------------------------------------------------------------------ +C Executable statements +C ------------------------------------------------------------------ + IF (IERROR.NE.0) THEN +C ------------------------------------------------------------------ +C Abnormal exit from calling routine +C ------------------------------------------------------------------ + IF (IFAIL.EQ.-1 .OR. IFAIL.EQ.0 .OR. IFAIL.EQ.-13 .OR. + & (IFAIL.GT.0 .AND. MOD(IFAIL/10,10).NE.0)) THEN +C ------------------------------------------------------------------ +C Noisy exit +C ------------------------------------------------------------------ + CALL X04AAF(0,NERR) + DO 20 I = 1, NREC + CALL X04BAF(NERR,REC(I)) +20 CONTINUE + IF (IFAIL.NE.-13) THEN + WRITE (MESS,FMT=99999) SRNAME, IERROR + CALL X04BAF(NERR,MESS) + IF (ABS(MOD(IFAIL,10)).NE.1) THEN +C ------------------------------------------------------------------ +C Hard failure +C ------------------------------------------------------------------ + CALL X04BAF(NERR, + & ' ** NAG HARD FAILURE - EXECUTION TERMINATED' + & ) + CALL P01ABZ + ELSE +C ------------------------------------------------------------------ +C Soft failure +C ------------------------------------------------------------------ + CALL X04BAF(NERR, + & ' ** NAG SOFT FAILURE - CONTROL RETURNED') + END IF + END IF + END IF + END IF + P01ABF = IERROR + RETURN + +99999 FORMAT (' ** ABNORMAL EXIT FROM NAG LIBRARY ROUTINE ',A,': IFAIL', + & ' =',I6) + + END +C ------------------------------------------------------------------ + SUBROUTINE P01ABZ +C ------------------------------------------------------------------ +C MARK 11.5 (F77) RELEASE. NAG COPYRIGHT 1986. +C +C Terminates execution when a hard failure occurs. +C +C ******************** IMPLEMENTATION NOTE ******************** +C THE FOLLOWING STOP STATEMENT MAY BE REPLACED BY A CALL TO AN +C IMPLEMENTATION-DEPENDENT ROUTINE TO DISPLAY A MESSAGE AND/OR +C TO ABORT THE PROGRAM. +C ************************************************************* +C ------------------------------------------------------------------ +C Executable statements +C ------------------------------------------------------------------ + STOP + END +C ------------------------------------------------------------------ + DOUBLE PRECISION FUNCTION X02AJF() +C ------------------------------------------------------------------ +C MARK 12 RELEASE. NAG COPYRIGHT 1986. +C +C Returns (1/2)*B**(1-P) if ROUNDS is .TRUE. +C Returns B**(1-P) otherwise +C ------------------------------------------------------------------ +C Executable statements +C ------------------------------------------------------------------ + X02AJF = 0.11102230246251568000D-015 + + RETURN + END +C ------------------------------------------------------------------ + DOUBLE PRECISION FUNCTION X02AKF() +C ------------------------------------------------------------------ +C MARK 12 RELEASE. NAG COPYRIGHT 1986. +C +C Returns B**(EMIN-1) (The smallest positive model number) +C ------------------------------------------------------------------ +C Executable statements +C ------------------------------------------------------------------ + X02AKF = 0.22250738585072014000D-307 + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE X04AAF(I,NERR) +C ------------------------------------------------------------------ +C MARK 7 RELEASE. NAG COPYRIGHT 1978. +C MARK 7C REVISED IER-190 (MAY 1979). +C MARK 11.5 (F77) REVISED (SEPT 1985). +C +C If I = 0, sets NERR to current error message unit number +C (stored in NERR1). +C If I = 1, changes current error message unit number to +C value specified by NERR. +C ------------------------------------------------------------------ +C Scalar arguments +C ------------------------------------------------------------------ + INTEGER I, NERR +C ------------------------------------------------------------------ +C Local scalars +C ------------------------------------------------------------------ + INTEGER NERR1 +C ------------------------------------------------------------------ +C Save statement +C ------------------------------------------------------------------ + SAVE NERR1 +C ------------------------------------------------------------------ +C Data statements +C ------------------------------------------------------------------ + DATA NERR1/0/ +C ------------------------------------------------------------------ +C Executable statements +C ------------------------------------------------------------------ + IF (I.EQ.0) NERR = NERR1 + IF (I.EQ.1) NERR1 = NERR + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE X04BAF(NOUT,REC) +C ------------------------------------------------------------------ +C MARK 11.5 (F77) RELEASE. NAG COPYRIGHT 1986. +C +C X04BAF writes teh contents of REC to the unit dfeined by NOUT. +C +C Trailing blanks are not output, except that if REC is entirely +C blank, a single blank character is output. +C If NOUT.LT.0, i.e. if NOUT is not a valid FORTRAN unit identifier, +C then no output occurs. +C ------------------------------------------------------------------ +C Scalar arguments +C ------------------------------------------------------------------ + INTEGER NOUT + CHARACTER*(*) REC +C ------------------------------------------------------------------ +C Local scalars +C ------------------------------------------------------------------ + INTEGER I +C ------------------------------------------------------------------ +C Intrinsic functions +C ------------------------------------------------------------------ + INTRINSIC LEN +C ------------------------------------------------------------------ +C Executable statements +C ------------------------------------------------------------------ + IF (NOUT.GE.0) THEN +C ------------------------------------------------------------------ +C Remove trailing blanks +C ------------------------------------------------------------------ + DO 20 I = LEN(REC), 2, -1 + IF (REC(I:I).NE.' ') GO TO 40 +20 CONTINUE +C ------------------------------------------------------------------ +C Write record to external file +C ------------------------------------------------------------------ +40 WRITE (NOUT,FMT=99999) REC(1:I) + END IF + RETURN + +99999 FORMAT (A) + + END diff -u -N -r xplor-nih-2.9.2/source/fantaxplor.fcm PARAxplor-nih-2.9.2/source/fantaxplor.fcm --- xplor-nih-2.9.2/source/fantaxplor.fcm 1970-01-01 01:00:00.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/fantaxplor.fcm 2003-12-05 14:34:36.000000000 +0100 @@ -0,0 +1,8 @@ +C ------------------------------------------------------------------ +C FANTALIN additional stuff +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + CHARACTER*4 NRESFA,NTYPEFA,NRESIDFA + + COMMON NRESFA(1000),NTYPEFA(1000),NRESIDFA(1000) diff -u -N -r xplor-nih-2.9.2/source/sup2.fcm PARAxplor-nih-2.9.2/source/sup2.fcm --- xplor-nih-2.9.2/source/sup2.fcm 1970-01-01 01:00:00.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/sup2.fcm 2003-12-05 14:34:36.000000000 +0100 @@ -0,0 +1,72 @@ + IMPLICIT NONE + INTEGER MAXOS,MAXFE,MAXSTR,MAXDEF,MPAR,MP,NP,MH + + PARAMETER (MAXOS=1000, MAXFE=2, MAXSTR=2, MAXDEF=1000, + * MPAR=5, MP=MAXFE*MPAR+1, NP=MAXFE*MPAR, MH=MAXOS*MAXSTR) + + CHARACTER FILENAME3*132 + + + DOUBLE PRECISION CX,CY,CZ, + * FX,FY,FZ, + * RCCPTX,RCCPTY,RCCPTZ, + * CENTER,SSX,SSY,SSZ, + * TSX,TSY,TSZ + + DOUBLE PRECISION TOLPROT,WPROT,SHIFT, + * OBS,MLPROT, + * PHI,TETA,OMEGA, + * A1DIP,A2DIP, + * OPTPHI,OPTTETA,OPTOMEGA, + * OPTA1,OPTA2,STR1,STR2, + * STR3,STR4,STR5,STR6, + * STR7,STR8, + * ASXX, ASYY, ASZZ, + * AXSYS, AYSYS, AZSYS, + * TOLDIP,RESID,OLDRESID,PERC,X1O,Y1O,Z1O,X2O + * ,Y2O,Z2O,X3O,Y3O,Z3O,DELTA,SIMP,RAND,EXTR + + DOUBLE PRECISION VAL,CASH,ZRR,X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,X4,Y4 + + DOUBLE PRECISION Z4,TETA2,PHI2,P,T,O,R,AL,BE,GA,RNDOM,PERM,Y + DOUBLE PRECISION APP,AVERAGE,ASTDEV,ERR,PBAR,PR,YR,PK,YK,PE,YE + + DOUBLE PRECISION PT,YT,PC,YC,OLDP,OLDPOPPA,RMR,TMP1,VETT,A1D,A2D + DOUBLE PRECISION TMP2,COSTMULT,PRIMO,TETA1,PHI1,G1,G2,G3,G4,G5,RM + DOUBLE PRECISION SMED,STDEV,VV,VMAX,VMIN,PI,A1,A2,SAXSYS,SAYSYS, + * SAZSYS,SAXX,SAYY,SAZZ,SDZX,APPX,APPY,SXX,SYY,SXZ,SYX,SXY, + * SYZ,SZX,SZY,SZZ,VZX,VZY,VZZ,CDIR,CHIPARA,CHIPERP + + + INTEGER IOLDVIO,IVIOLATION,INTSYS, + * NHP,NFE,NSTR,NUMS,NUMD,NUMT,NAT,NSYSTEM,NGRID + * ,IHP,NIA,L,NK,I,J,ITER,IL,IM,JM,N,M,NPROTML + * ,MED,NORES,NUMBERATOM,NR,MDIPO + + + DIMENSION XP(MAXSTR,MAXDEF),YP(MAXSTR,MAXDEF),ZP(MAXSTR,MAXDEF) + DIMENSION NUM_RES(MAXSTR,MAXDEF) + + DIMENSION NUMRES(MH),R(MAXSTR,MAXOS),TETA1(MAXSTR,MAXOS), + * PHI1(MAXSTR,MAXOS) + + COMMON /SIMPLEX/ + * CX(MAXSTR,MAXOS),CY(MAXSTR,MAXOS),CZ(MAXSTR,MAXOS), + * FX(MAXSTR,MAXOS),FY(MAXSTR,MAXOS),FZ(MAXSTR,MAXOS), + * RCCPTX(MAXFE,3), RCCPTY(MAXFE,3),RCCPTZ(MAXFE,3), + * CENTER(MAXFE,3),SSX(MAXFE,3),SSY(MAXFE,3),SSZ(MAXFE,3), + * TSX(MAXFE,3),TSY(MAXFE,3),TSZ(MAXFE,3), + * TOLPROT(MH),WPROT(MH),SHIFT(MH), + * OBS(MH),MLPROT(MH), + * PHI(MAXFE),TETA(MAXFE),OMEGA(MAXFE), + * A1DIP(MAXFE),A2DIP(MAXFE), + * OPTPHI(MAXFE),OPTTETA(MAXFE),OPTOMEGA(MAXFE), + * OPTA1(MAXFE),OPTA2(MAXFE),STR1(MAXSTR),STR2(MAXSTR), + * STR3(MAXSTR),STR4(MAXSTR),STR5(MAXSTR),STR6(MAXSTR), + * STR7(MAXSTR),STR8(MAXSTR), + * ASXX(MAXFE,3), ASYY(MAXFE,3), ASZZ(MAXFE,3), + * AXSYS(MAXFE,3), AYSYS(MAXFE,3), AZSYS(MAXFE,3), + * TOLDIP,RESID,OLDRESID,IOLDVIO,IVIOLATION,INTSYS, + * NHP,NFE,NSTR,NUMS,NUMD,NUMT,NAT,NSYSTEM,NGRID, + * FILENAME3,MDIPO,CHIPARA,CHIPERP + diff -u -N -r xplor-nih-2.9.2/source/supccr.fcm PARAxplor-nih-2.9.2/source/supccr.fcm --- xplor-nih-2.9.2/source/supccr.fcm 1970-01-01 01:00:00.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/supccr.fcm 2003-12-05 14:34:36.000000000 +0100 @@ -0,0 +1,72 @@ + IMPLICIT NONE + INTEGER MAXOS,MAXFE,MAXSTR,MAXDEF,MPAR,MP,NP,MH + + PARAMETER (MAXOS=1000, MAXFE=2, MAXSTR=2, MAXDEF=1000, + * MPAR=1, MP=MAXFE*MPAR+1, NP=MAXFE*MPAR, MH=MAXOS*MAXSTR) + + CHARACTER FILENAME3*132 + + + DOUBLE PRECISION CX,CY,CZ, + * FX,FY,FZ, + * RCCPTX,RCCPTY,RCCPTZ, + * CENTER,SSX,SSY,SSZ, + * TSX,TSY,TSZ + + DOUBLE PRECISION TOLPROT,WPROT,CSHIFT, + * OBS,MLPROT, + * PHI,TETA,OMEGA, + * A1DIP,A2DIP, + * OPTPHI,OPTTETA,OPTOMEGA, + * OPTA1,OPTA2,STR1,STR2, + * STR3,STR4,STR5,STR6, + * STR7,STR8, + * ASXX, ASYY, ASZZ, + * AXSYS, AYSYS, AZSYS, + * TOLDIP,RESID,OLDRESID,PERC,X1O,Y1O,Z1O,X2O + * ,Y2O,Z2O,X3O,Y3O,Z3O,DELTA,SIMP,RAND,EXTR + DOUBLE PRECISION AK1,RK1,RK2 + DOUBLE PRECISION VAL,CCASH,ZRR,X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,X4,Y4 + + DOUBLE PRECISION Z4,TETA2,PHI2,P,T,O,R,AL,BE,GA,RNDOM,PERM,Y + DOUBLE PRECISION APP,AVERAGE,ASTDEV,ERR,PBAR,PR,YR,PK,YK,PE,YE + + DOUBLE PRECISION PT,YT,PC,YC,OLDP,OLDPOPPA,RMR,TMP1,VETT,A1D,A2D + DOUBLE PRECISION TMP2,COSTMULT,PRIMO,TETA1,PHI1,G1,G2,G3,G4,G5,RM + DOUBLE PRECISION SMED,STDEV,VV,VMAX,VMIN,PI,A1,A2,SAXSYS,SAYSYS, + * SAZSYS,SAXX,SAYY,SAZZ,SDZX,APPX,APPY,SXX,SYY,SXZ,SYX,SXY, + * SYZ,SZX,SZY,SZZ,VZX,VZY,VZZ,CDIR + + + INTEGER IOLDVIO,IVIOLATION,INTSYS, + * NHP,NFE,NSTR,NUMS,NUMD,NUMT,NAT,NSYSTEM,NGRID + * ,IHP,NIA,L,NK,I,J,ITER,IL,IM,JM,N,M,NPROTML + * ,MED,NORES,NUMBERATOM,NR,MDIPO + + + DIMENSION XP(MAXSTR,MAXDEF),YP(MAXSTR,MAXDEF),ZP(MAXSTR,MAXDEF) + DIMENSION NUM_RES(MAXSTR,MAXDEF) + + DIMENSION NUMRES(MH),R(MAXSTR,MAXOS),TETA1(MAXSTR,MAXOS), + * PHI1(MAXSTR,MAXOS) + + COMMON /SIMPLEX/ + * CX(MAXSTR,MAXOS),CY(MAXSTR,MAXOS),CZ(MAXSTR,MAXOS), + * FX(MAXSTR,MAXOS),FY(MAXSTR,MAXOS),FZ(MAXSTR,MAXOS), + * RCCPTX(MAXFE,3), RCCPTY(MAXFE,3),RCCPTZ(MAXFE,3), + * CENTER(MAXFE,3),SSX(MAXFE,3),SSY(MAXFE,3),SSZ(MAXFE,3), + * TSX(MAXFE,3),TSY(MAXFE,3),TSZ(MAXFE,3), + * TOLPROT(MH),WPROT(MH),CSHIFT(MH), + * OBS(MH),MLPROT(MH),AK1(MAXOS),RK1(MAXOS),RK2(MAXOS), + * PHI(MAXFE),TETA(MAXFE),OMEGA(MAXFE), + * A1DIP(MAXFE),A2DIP(MAXFE), + * OPTPHI(MAXFE),OPTTETA(MAXFE),OPTOMEGA(MAXFE), + * OPTA1(MAXFE),OPTA2(MAXFE),STR1(MAXSTR),STR2(MAXSTR), + * STR3(MAXSTR),STR4(MAXSTR),STR5(MAXSTR),STR6(MAXSTR), + * STR7(MAXSTR),STR8(MAXSTR), + * ASXX(MAXFE,3), ASYY(MAXFE,3), ASZZ(MAXFE,3), + * AXSYS(MAXFE,3), AYSYS(MAXFE,3), AZSYS(MAXFE,3), + * TOLDIP,RESID,OLDRESID,IOLDVIO,IVIOLATION,INTSYS, + * NHP,NFE,NSTR,NUMS,NUMD,NUMT,NAT,NSYSTEM,NGRID, + * FILENAME3,MDIPO + diff -u -N -r xplor-nih-2.9.2/source/xangle.f PARAxplor-nih-2.9.2/source/xangle.f --- xplor-nih-2.9.2/source/xangle.f 1970-01-01 01:00:00.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/xangle.f 2003-12-05 14:34:36.000000000 +0100 @@ -0,0 +1,772 @@ +C ------------------------------------------------------------------ + SUBROUTINE EXANGLES (EANGLES, WHICH) +C ------------------------------------------------------------------ +C Calls EXANGLES2, which does the actual energy calculation +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "xangle.fcm" + INCLUDE "heap.fcm" + + DOUBLE PRECISION EANGLES + CHARACTER*7 WHICH + + CALL EXANGLES2(EANGLES,HEAP(XANGLESJPTR),HEAP(XANGLESKPTR), + & HEAP(XANGLESJLST),HEAP(XANGLESKLST), + & HEAP(XANGLESOBS1PTR),HEAP(XANGLESOBS2PTR), + & HEAP(XANGLESERRPTR), + & HEAP(CALCXANGLESPTR),WHICH) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE EXANGLES2 (EANGLES,ATOMJPTR,ATOMKPTR,ATOMJLST,ATOMKLST, + & XANGLESOBS1,XANGLESOBS2,XANGLESERR, + & XANGLESCALC,WHICH) +C ------------------------------------------------------------------ +C Calculates RDC-ANGLES energies +C +C Energies are of the form: +C E = K1*(SIN(ANG)**2) +C where: +C K1 = Force constant, +C ANG = The angle between the target vector and the actual vector +C +C Note that: +C Atom J = Detected nucleus (e.g. H) +C Atom K = Coupled nucleus (e.g. N) +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "psf.fcm" + INCLUDE "coord.fcm" + INCLUDE "numbers.fcm" + INCLUDE "deriv.fcm" + INCLUDE "xangle.fcm" + INCLUDE "consta.fcm" + + INTEGER ATOMJPTR(*),ATOMKPTR(*) + INTEGER ATOMJLST(*),ATOMKLST(*) + DOUBLE PRECISION XANGLESOBS1(*),XANGLESOBS2(*),XANGLESERR(*) + DOUBLE PRECISION XANGLESCALC(*),EANGLES + CHARACTER*7 WHICH +C ------------------------------------------------------------------ +C Local variables +C ------------------------------------------------------------------ + INTEGER COUNT,CLASS,MM,NN,II,I,J,K,L,NA,M,N + PARAMETER (NA=20) + INTEGER DMM,DNN,OUTFLAG(NA),MCOUNT,NCOUNT,PCOUNT + DOUBLE PRECISION XJ,XK,YJ,YK,ZJ,ZK,CALCXANGLES,E, + & DFXJ(NA),DFYJ(NA),DFZJ(NA), + & DFXK(NA),DFYK(NA),DFZK(NA),DF1, + & O1,O2,O3,O4,O5,O6,O7,O8,O9,O10,O11,O12, + & PS,DJK,DJK2,DJK3 + DOUBLE PRECISION OBSXANGLES1,OBSXANGLES2,ERRXANGLES,K1, + & A,B, + & DT,DP,DELTA,DELTAXANGLES, + & DD +C ------------------------------------------------------------------ +C Zero out partial energy +C ------------------------------------------------------------------ + EANGLES = ZERO + + CLASS = 1 + K1=XANGLESFORCES(1) + + COUNT=0 + + DO WHILE(COUNT.LT.ANGLESNUM) + COUNT = COUNT+1 +C ------------------------------------------------------------------ +C Reset individual E to zero +C ------------------------------------------------------------------ + E=0 + + DO WHILE ((XANGLESASSNDX(CLASS).LT.COUNT).OR. + & (XANGLESASSNDX(CLASS).EQ.0)) + CLASS=CLASS+1 + END DO + + IF (XANGLESASSNDX(CLASS).EQ.XANGLESASSNDX(CLASS-1)) THEN + COUNT=COUNT-1 + ENDIF + + K1=XANGLESFORCES(CLASS) +C ------------------------------------------------------------------ +C Note there should only be one atom for J and K +C ------------------------------------------------------------------ + J=ATOMJPTR(COUNT)+1 + K=ATOMKPTR(COUNT)+1 + + XJ=X(ATOMJLST(J)) + XK=X(ATOMKLST(K)) + YJ=Y(ATOMJLST(J)) + YK=Y(ATOMKLST(K)) + ZJ=Z(ATOMJLST(J)) + ZK=Z(ATOMKLST(K)) + + OBSXANGLES1=XANGLESOBS1(COUNT) + OBSXANGLES2=XANGLESOBS2(COUNT) + ERRXANGLES=XANGLESERR(COUNT) +C ------------------------------------------------------------------ +C Initialize calculated RDC-angles and counter +C ------------------------------------------------------------------ + CALCXANGLES=0 + II=0 + MCOUNT=0 + NCOUNT=0 +C ------------------------------------------------------------------ +C Check for correct permutations of paired atoms to get the nucleus- +C nucleus vector. This depends solely on the order of the assigned +C atoms. OUTFLAG=1 indicates that the permutation is not allowed. +C ------------------------------------------------------------------ + DMM=ATOMJPTR(COUNT+1)-ATOMJPTR(COUNT) + DNN=ATOMKPTR(COUNT+1)-ATOMKPTR(COUNT) + DO MM=ATOMJPTR(COUNT)+1,ATOMJPTR(COUNT+1) + MCOUNT=MCOUNT+1 + NCOUNT=0 + DO NN=ATOMKPTR(COUNT)+1,ATOMKPTR(COUNT+1) + NCOUNT=NCOUNT+1 + II=II+1 + IF ((DMM.GT.1).AND.(DNN.GT.1)) THEN + IF (MCOUNT.EQ.NCOUNT) THEN + OUTFLAG(II)=0 + ELSE + OUTFLAG(II)=1 + ENDIF + ELSE + OUTFLAG(II)=0 + END IF + END DO + END DO + + II=0 + PCOUNT=0 + DO MM=ATOMJPTR(COUNT)+1,ATOMJPTR(COUNT+1) + DO NN=ATOMKPTR(COUNT)+1,ATOMKPTR(COUNT+1) + II=II+1 + IF (OUTFLAG(II).NE.1) THEN + PCOUNT=PCOUNT+1 +C ------------------------------------------------------------------ +C Calculate the function and its derivatives +C ------------------------------------------------------------------ + OBSXANGLES1=OBSXANGLES1*PI/180. + IF (OBSXANGLES1.LT.0.) OBSXANGLES1=OBSXANGLES1+2.*PI + OBSXANGLES2=OBSXANGLES2*PI/180. + IF (OBSXANGLES2.LT.0.) OBSXANGLES2=OBSXANGLES2+2.*PI + + O1=SIN(OBSXANGLES1)*COS(OBSXANGLES2) + O2=SIN(OBSXANGLES1)*SIN(OBSXANGLES2) + O3=COS(OBSXANGLES1) + O4=XJ-XK + O5=YJ-YK + O6=ZJ-ZK + O7=O4**2 + O8=O5**2 + O9=O6**2 + DJK2=O7+O8+O9 + DJK=SQRT(DJK2) + DJK3=DJK**3 + PS=O1*O4+O2*O5+O3*O6 + O10=O1*DJK-PS*O4/DJK + O11=O2*DJK-PS*O5/DJK + O12=O3*DJK-PS*O6/DJK + + DFXJ(II)=-2.*PS*O10/DJK3 + DFYJ(II)=-2.*PS*O11/DJK3 + DFZJ(II)=-2.*PS*O12/DJK3 + DFXK(II)=-DFXJ(II) + DFYK(II)=-DFYJ(II) + DFZK(II)=-DFZJ(II) +C ------------------------------------------------------------------ +C Calculate the cosine square of the angle +C ------------------------------------------------------------------ + CALCXANGLES=(PS/DJK)**2 + END IF + END DO + END DO + + IF (WHICH.EQ.'ANALYZE') THEN + XANGLESCALC(COUNT)=CALCXANGLES + END IF +C ------------------------------------------------------------------ +C The deviation is the sine square of the angle +C ------------------------------------------------------------------ + DELTA=(1.-CALCXANGLES) +C ------------------------------------------------------------------ +C Adjust the deviation based on the error range (expressed as the +C accepted deviation of the cosine square from 1) +C ------------------------------------------------------------------ + IF ((DELTA.GT.0.000).AND.(DELTA.GT.ERRXANGLES)) THEN + DELTAXANGLES=DELTA-ERRXANGLES + ELSE + DELTAXANGLES=0.0 + END IF + + DD=DELTAXANGLES + + II=0 + PCOUNT=0 + DO MM=ATOMJPTR(COUNT)+1,ATOMJPTR(COUNT+1) + DO NN=ATOMKPTR(COUNT)+1,ATOMKPTR(COUNT+1) + II=II+1 + IF (OUTFLAG(II).NE.1) THEN + PCOUNT=PCOUNT+1 + IF (DD.EQ.0.0) THEN + E=E+0.0 + DF1=0.0 + ELSE + E=E+K1*DELTAXANGLES + DF1=K1 + END IF + END IF + END DO + END DO +C ------------------------------------------------------------------ +C Accumulate energy +C ------------------------------------------------------------------ + EANGLES=EANGLES+E +C ------------------------------------------------------------------ +C Now update forces if in energy/force mode +C ------------------------------------------------------------------ + IF (WHICH.NE.'ANALYZE') THEN + II=0 + DO MM=ATOMJPTR(COUNT)+1,ATOMJPTR(COUNT+1) + DO NN=ATOMKPTR(COUNT)+1,ATOMKPTR(COUNT+1) + II=II+1 + IF (OUTFLAG(II).NE.1) THEN + DX(ATOMJLST(J))=DX(ATOMJLST(J))+DF1*DFXJ(II) + DY(ATOMJLST(J))=DY(ATOMJLST(J))+DF1*DFYJ(II) + DZ(ATOMJLST(J))=DZ(ATOMJLST(J))+DF1*DFZJ(II) + DX(ATOMKLST(K))=DX(ATOMKLST(K))+DF1*DFXK(II) + DY(ATOMKLST(K))=DY(ATOMKLST(K))+DF1*DFYK(II) + DZ(ATOMKLST(K))=DZ(ATOMKLST(K))+DF1*DFZK(II) + END IF + END DO + END DO + END IF + END DO + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE READXANGLES +C ------------------------------------------------------------------ +C Reads in RDC-ANGLES information +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "comand.fcm" + INCLUDE "xangle.fcm" + INCLUDE "funct.fcm" + INCLUDE "psf.fcm" + INCLUDE "heap.fcm" + INCLUDE "numbers.fcm" + + INTEGER COUNT,SPTR,OLDCLASS,OLDMAXXANGLES + DOUBLE PRECISION K1,CUTOFF + CHARACTER*4 THENAME + + SPTR=ALLHP(INTEG4(NATOM)) + CALL PUSEND('XANGLE>') +862 CONTINUE + CALL NEXTWD('XANGLE>') + CALL MISCOM('XANGLE>',USED) + IF (.NOT.USED) THEN +C ------------------------------------------------------------------ +C Documentation +C ------------------------------------------------------------------ + IF (WD(1:4).EQ.'HELP') THEN + WRITE(DUNIT,'(10X,A)') + & ' XANGLE {} END ', + & ' :== ', + & ' ASSIGN ', + & ' * Restraint: Detected_nucleus Coupled_nucleus', + & ' Theta Phi Err *', + & ' CLASsification ', + & ' * Starts a new class *', + & ' FORCeconstant ', + & ' * Force constant for the current class *', + & ' NREStraints ', + & ' * Number of slots for restraints to allocate *', + & ' PRINt THREshold ', + & ' * Prints violations larger than the THREshold *', + & ' RESEt', + & ' * Erases the restraint table, but keeps NRES *' +C ------------------------------------------------------------------ +C Get class name. Determine if it's an already-defined class. +C Insert a new class if it's not. +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'CLAS') THEN + OLDCLASS=CURCLASS + CALL NEXTA4('Class name =',THENAME) + MODE=NEW + DO COUNT=1,NCLASSES + IF (XANGLESCLASSNAMES(COUNT).EQ.THENAME) THEN + MODE=UPDATE + CURCLASS=COUNT + END IF + END DO + IF (MODE.EQ.NEW) THEN +C ------------------------------------------------------------------ +C Make sure you can't add more than the maximum number of classes +C ------------------------------------------------------------------ + IF (OLDCLASS.EQ.MAXXANGLESCLASSES) THEN + CALL DSPERR('XANGLE','Too many classes.') + CALL DSPERR('XANGLE', + & 'Increase MAXXANGLESCLASSES and recompile.') + CALL WRNDIE(-5, 'READXANGLES', + & 'Too many RDC-angles classes.') + END IF + NCLASSES=NCLASSES+1 + CURCLASS=NCLASSES + XANGLESCLASSNAMES(CURCLASS)=THENAME +C ------------------------------------------------------------------ +C If this isn't the first class, close off the old class +C ------------------------------------------------------------------ + IF (NCLASSES.GT.1) THEN + XANGLESASSNDX(OLDCLASS)=ANGLESNUM + END IF + END IF +C ------------------------------------------------------------------ +C Set force constant for current class +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'FORC') THEN +C ------------------------------------------------------------------ +C Start a default class if there isn't one defined +C ------------------------------------------------------------------ + IF (CURCLASS.EQ.0) THEN + NCLASSES = 1 + CURCLASS = 1 + END IF + CALL NEXTF('Force constant =',K1) + WRITE(DUNIT,'(A,A,A,F8.3)') + & 'Setting force constant for class ', + & XANGLESCLASSNAMES(CURCLASS),' to ',K1 + XANGLESFORCES(CURCLASS)=K1 +C ------------------------------------------------------------------ +C Reset RDC-angles database +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'RESE') THEN + CALL XANGLESDEFAULTS + CALL ALLOCXANGLES(0,MAXXANGLES) +C ------------------------------------------------------------------ +C Change number of assignment slots +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'NRES') THEN + OLDMAXXANGLES=MAXXANGLES + CALL NEXTI('Number of slots =',MAXXANGLES) + CALL ALLOCXANGLES(OLDMAXXANGLES,MAXXANGLES) + WRITE(DUNIT,'(A,I8,A)') + & 'XANGLE: Allocating space for',MAXXANGLES, + & 'number of restraints.' +C ------------------------------------------------------------------ +C Read in an assignment +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'ASSI') THEN +C ------------------------------------------------------------------ +C Make sure you can't add more assignments than you have slots for +C ------------------------------------------------------------------ + IF (XANGLESMX.EQ.MAXXANGLES) THEN + CALL DSPERR('XANGLE','Too many assignments.') + CALL DSPERR('XANGLE', + & 'Increasing NREStraints by 100.') + OLDMAXXANGLES=MAXXANGLES + MAXXANGLES=MAXXANGLES+100 + CALL ALLOCXANGLES(OLDMAXXANGLES,MAXXANGLES) + END IF +C ------------------------------------------------------------------ +C If there isn't a class specified, start a default class +C ------------------------------------------------------------------ + IF (CURCLASS.EQ.0) THEN + NCLASSES=1 + CURCLASS=1 + END IF + CALL READXANGLES2(HEAP(XANGLESJPTR),HEAP(XANGLESKPTR), + & HEAP(XANGLESJLST),HEAP(XANGLESKLST), + & HEAP(XANGLESOBS1PTR), + & HEAP(XANGLESOBS2PTR), + & HEAP(XANGLESERRPTR),HEAP(SPTR)) +C ------------------------------------------------------------------ +C Print violations +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'PRIN') THEN + CALL NEXTWD('PRINT>') + IF (WD(1:4).NE.'THRE') THEN + CALL DSPERR('XANGLE', + & 'PRINt expects THREshold parameter.') + ELSE + CALL NEXTF('THREshold =',CUTOFF) + IF (CUTOFF.LT.ZERO) THEN + CALL DSPERR('XANGLE', + & 'Cutoff must be positive.') + CUTOFF = ABS(CUTOFF) + END IF + CALL NEXTA4('ALL OR CLASs>',THENAME) + IF (THENAME(1:3).EQ.'ALL') THEN + DO COUNT=1,NCLASSES + PRINTCLASS(COUNT)=.TRUE. + END DO + ELSE IF (THENAME(1:4).EQ.'CLAS') THEN + CALL NEXTA4('CLASS NAME =',THENAME) + DO COUNT=1,NCLASSES + IF (XANGLESCLASSNAMES(COUNT).EQ.THENAME) THEN + PRINTCLASS(COUNT)=.TRUE. + ELSE + PRINTCLASS(COUNT)=.FALSE. + END IF + END DO + ELSE + DO COUNT=1,NCLASSES + PRINTCLASS(COUNT)=.TRUE. + END DO + END IF + CALL PRINTXANGLES(CUTOFF,HEAP(CALCXANGLESPTR), + & HEAP(XANGLESOBS1PTR), + & HEAP(XANGLESOBS2PTR), + & HEAP(XANGLESERRPTR), + & HEAP(XANGLESJPTR), + & HEAP(XANGLESKPTR), + & HEAP(XANGLESJLST), + & HEAP(XANGLESKLST)) + END IF +C ------------------------------------------------------------------ +C Check for END statement +C ------------------------------------------------------------------ + ELSE + CALL CHKEND('XANGLE>', DONE) + END IF + END IF + IF (.NOT.(DONE)) GOTO 862 + DONE=.FALSE. + CALL FREHP(SPTR,INTEG4(NATOM)) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE ALLOCXANGLES (OLDSIZE, NEWSIZE) +C ------------------------------------------------------------------ +C Resets RDC-angles arrays to hold size entries +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "funct.fcm" + INCLUDE "xangle.fcm" + + INTEGER OLDSIZE,NEWSIZE + + IF (OLDSIZE.NE.0) THEN + CALL FREHP(XANGLESJPTR,INTEG4(OLDSIZE)) + CALL FREHP(XANGLESKPTR,INTEG4(OLDSIZE)) + + CALL FREHP(XANGLESJLST,INTEG4(OLDSIZE)) + CALL FREHP(XANGLESKLST,INTEG4(OLDSIZE)) + + CALL FREHP(XANGLESOBS1PTR,IREAL8(OLDSIZE)) + CALL FREHP(XANGLESOBS2PTR,IREAL8(OLDSIZE)) + CALL FREHP(XANGLESERRPTR,IREAL8(OLDSIZE)) + CALL FREHP(CALCXANGLESPTR,IREAL8(OLDSIZE)) + END IF + + XANGLESJPTR=ALLHP(INTEG4(NEWSIZE)) + XANGLESKPTR=ALLHP(INTEG4(NEWSIZE)) + + XANGLESJLST=ALLHP(INTEG4(NEWSIZE)) + XANGLESKLST=ALLHP(INTEG4(NEWSIZE)) + + XANGLESOBS1PTR=ALLHP(IREAL8(NEWSIZE)) + XANGLESOBS2PTR=ALLHP(IREAL8(NEWSIZE)) + XANGLESERRPTR=ALLHP(IREAL8(NEWSIZE)) + CALCXANGLESPTR=ALLHP(IREAL8(NEWSIZE)) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE XANGLESDEFAULTS +C ------------------------------------------------------------------ +C Sets up defaults +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "xangle.fcm" + + INTEGER COUNT + + MODE=NEW + MAXXANGLES=200 + XANGLESMX=200 + ANGLESNUM=0 + NCLASSES=0 + CURCLASS=0 + DO COUNT=1,MAXXANGLESCLASSES + XANGLESCLASSNAMES(COUNT)='DEFAULT' + XANGLESASSNDX(COUNT)=0 + XANGLESFORCES(COUNT)=5.0 + END DO + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE READXANGLES2(ATOMJPTR,ATOMKPTR,ATOMJLST,ATOMKLST, + & XANGLESOBS1,XANGLESOBS2,XANGLESERR,SEL) +C ------------------------------------------------------------------ +C Reads actual RDC-angles assignments into arrays +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "psf.fcm" + INCLUDE "coord.fcm" + INCLUDE "xangle.fcm" + + INTEGER ATOMJPTR(*),ATOMKPTR(*),ATOMJLST(*),ATOMKLST(*),SEL(*) + INTEGER JTMP,KTMP,LTMP,II + DOUBLE PRECISION XANGLESOBS1(*),XANGLESOBS2(*),XANGLESERR(*) +C ------------------------------------------------------------------ +C Local variables +C ------------------------------------------------------------------ + INTEGER NSEL,INSERTPOS,COUNT,CURSTOP,OTHERSTOP,NFLAG + INTEGER I + DOUBLE PRECISION XANGLESO1,XANGLESO2,XANGLESE +C ------------------------------------------------------------------ +C If we're in UPDATE mode, make a space for the new line +C ------------------------------------------------------------------ + NFLAG = 0 + IF (MODE.EQ.UPDATE) THEN + DO COUNT=ANGLESNUM+1,XANGLESASSNDX(CURCLASS)+1,-1 + ATOMJPTR(COUNT)=ATOMJPTR(COUNT-1) + ATOMKPTR(COUNT)=ATOMKPTR(COUNT-1) + XANGLESOBS1(COUNT)=XANGLESOBS1(COUNT-1) + XANGLESOBS2(COUNT)=XANGLESOBS2(COUNT-1) + XANGLESERR(COUNT)=XANGLESERR(COUNT-1) + END DO + CURSTOP=XANGLESASSNDX(CURCLASS) + DO COUNT=1,NCLASSES + OTHERSTOP=XANGLESASSNDX(COUNT) + IF (OTHERSTOP.GT.CURSTOP) THEN + XANGLESASSNDX(COUNT)=OTHERSTOP+1 + END IF + END DO + XANGLESASSNDX(CURCLASS)=CURSTOP+1 + INSERTPOS=CURSTOP + ANGLESNUM=ANGLESNUM + 1 + ELSE + ANGLESNUM=ANGLESNUM + 1 + INSERTPOS=ANGLESNUM + XANGLESASSNDX(CURCLASS)=INSERTPOS + END IF +C ------------------------------------------------------------------ +C Reading in the atom selection in the constraint table +C ------------------------------------------------------------------ + CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.) + IF (NSEL.GT.1) THEN + CALL DSPERR('XANGLE', + & 'More than 1 atom in SEL for atom J. Using first') + END IF + IF (NSEL.EQ.0) THEN + CALL DSPERR('XANGLE','Error with atom selection') + NFLAG=1 + END IF + CALL MAKIND(SEL, NATOM, NSEL) + IF (INSERTPOS.EQ.1) ATOMJPTR(INSERTPOS)=0 + II=ATOMJPTR(INSERTPOS) + II=II+1 + IF (II.GT.XANGLESMX) XANGLESMX=II + ATOMJLST(II)=SEL(1) + ATOMJPTR(INSERTPOS+1)=II + + CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.) + IF (NSEL.GT.1) THEN + CALL DSPERR('XANGLE', + & 'More than 1 atom in SEL for atom K. Using first') + END IF + IF (NSEL.EQ.0) THEN + CALL DSPERR('XANGLE','Error with atom selection') + NFLAG=1 + END IF + CALL MAKIND(SEL,NATOM,NSEL) + IF (INSERTPOS.EQ.1) ATOMKPTR(INSERTPOS)=0 + II=ATOMKPTR(INSERTPOS) + II=II+1 + IF (II.GT.XANGLESMX) XANGLESMX=II + ATOMKLST(II)=SEL(1) + ATOMKPTR(INSERTPOS+1)=II +C ------------------------------------------------------------------ +C Reading in the observed Theta and Phi +C ------------------------------------------------------------------ + CALL NEXTF('Observed THETA =',XANGLESO1) + CALL NEXTF('Observed PHI =',XANGLESO2) + CALL NEXTF('Error on sine square =',XANGLESE) + + XANGLESOBS1(INSERTPOS)=XANGLESO1 + XANGLESOBS2(INSERTPOS)=XANGLESO2 + XANGLESERR(INSERTPOS)=XANGLESE +C ------------------------------------------------------------------ +C Check for error in atom selection. If there is one then reset the +C counter for restraint +C ------------------------------------------------------------------ + IF (NFLAG.EQ.1) THEN + ANGLESNUM=ANGLESNUM-1 + END IF + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE XANGLESINIT +C ------------------------------------------------------------------ +C Initializes RDC-angles stuff +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "xangle.fcm" + + CALL XANGLESDEFAULTS + CALL ALLOCXANGLES(0,MAXXANGLES) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE PRINTXANGLES (CUTOFF,XANGLESCALC,XANGLESOBS1, + & XANGLESOBS2,XANGLESERR, + & ATOMJPTR,ATOMKPTR,ATOMJLST,ATOMKLST) +C ------------------------------------------------------------------ +C Prints RDC-angles with deviation of cosine square from 1 larger +C than cutoff, calculates RMSD and puts it into $RESULT +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "xangle.fcm" + INCLUDE "comand.fcm" + INCLUDE "psf.fcm" + INCLUDE "numbers.fcm" + + DOUBLE PRECISION CUTOFF,XANGLESCALC(*),XANGLESOBS1(*), + & XANGLESOBS2(*),XANGLESERR(*) + INTEGER ATOMJLST(*),ATOMKLST(*) + INTEGER ATOMJPTR(*),ATOMKPTR(*) +C ------------------------------------------------------------------ +C Local variables +C ------------------------------------------------------------------ + DOUBLE PRECISION CALCXANGLES,OBSXANGLES1,OBSXANGLES2, + & DELTAXANGLES,DELTA,DP + INTEGER COUNT,CLASS,J,K,NUM,II + DOUBLE PRECISION RMS,VIOLS,ERRXANGLES,XANGLESENERGY + DOUBLE COMPLEX DUMMY2 + LOGICAL PRINTTHISCLASS + + RMS=ZERO + VIOLS=ZERO + NUM=0 +C ------------------------------------------------------------------ +C Make sure that the array of calculated RDC-angles is up to date +C ------------------------------------------------------------------ + CALL EXANGLES(XANGLESENERGY,'ANALYZE') + WRITE (PUNIT,'(A)') 'The following RDC-angles have' + WRITE (PUNIT,'(A)') 'deviation of cosine square from 1' + WRITE (PUNIT,'(A)') 'larger than or equal to the cutoff:' +C ------------------------------------------------------------------ +C Write out first class heading +C ------------------------------------------------------------------ + CLASS=1 + PRINTTHISCLASS=PRINTCLASS(CLASS) + IF (PRINTTHISCLASS) THEN + WRITE (PUNIT,'(A,A)') 'Class ',XANGLESCLASSNAMES(1) + END IF +C ------------------------------------------------------------------ +C For every RDC-angles entry... +C ------------------------------------------------------------------ + COUNT=0 + DO WHILE(COUNT.LT.ANGLESNUM) + COUNT=COUNT+1 +C ------------------------------------------------------------------ +C Is this the start of a new class? +C ------------------------------------------------------------------ + IF (XANGLESASSNDX(CLASS).LT.COUNT) THEN + CLASS=CLASS+1 + PRINTTHISCLASS=PRINTCLASS(CLASS) + IF (XANGLESASSNDX(CLASS).EQ.XANGLESASSNDX(CLASS-1)) THEN + COUNT=COUNT-1 + END IF + IF (PRINTTHISCLASS) THEN + WRITE(PUNIT,'(A,A)') 'Class ',XANGLESCLASSNAMES(CLASS) + END IF + END IF +C ------------------------------------------------------------------ +C If this assignment is in a class to be printed +C and make sure there is an entry for that class +C ------------------------------------------------------------------ + IF((PRINTTHISCLASS).AND. + & (XANGLESASSNDX(CLASS).NE.XANGLESASSNDX(CLASS-1))) THEN +C ------------------------------------------------------------------ +C Always update RMSD +C ------------------------------------------------------------------ + CALCXANGLES=XANGLESCALC(COUNT) + OBSXANGLES1=XANGLESOBS1(COUNT) + OBSXANGLES2=XANGLESOBS2(COUNT) + ERRXANGLES=XANGLESERR(COUNT) + DP=(1.-CALCXANGLES) + + IF ((DP.GT.0.000).AND.(DP.GT.ERRXANGLES)) THEN + DELTAXANGLES=DP-ERRXANGLES + ELSE + DELTAXANGLES=0.0 + END IF + + RMS=RMS+DELTAXANGLES**2 + NUM=NUM+1 +C ------------------------------------------------------------------ +C Print out deviations larger than cutoff +C and update number of violations +C ------------------------------------------------------------------ + IF (ABS(DELTAXANGLES).GE.CUTOFF) THEN + J=ATOMJLST(ATOMJPTR(COUNT)+1) + K=ATOMKLST(ATOMKPTR(COUNT)+1) + WRITE(PUNIT,'(A,A)') '==============================', + & '==============================' + WRITE(PUNIT,'(A)') ' Set-J-atoms' + DO II=ATOMJPTR(COUNT)+1,ATOMJPTR(COUNT+1) + J=ATOMJLST(II) + WRITE(PUNIT,'(9X,4(1X,A))') SEGID(J),RESID(J), + & RES(J),TYPE(J) + END DO + WRITE(PUNIT,'(A)') ' Set-K-atoms' + DO II=ATOMKPTR(COUNT)+1,ATOMKPTR(COUNT+1) + K=ATOMKLST(II) + WRITE(PUNIT,'(9X,4(1X,A))') SEGID(K),RESID(K), + & RES(K),TYPE(K) + END DO + WRITE(PUNIT,'(2(2X,A,1X,F8.3))') + & 'Calc: ',CALCXANGLES,'Obs: ',1. + WRITE(PUNIT,'(2X,A,1X,F8.3,2X,A,1X,F8.3)') + & 'Error: ',ERRXANGLES,'Delta: ',DELTAXANGLES + VIOLS=VIOLS+ONE + END IF + END IF + END DO + + IF (NUM.GT.0) THEN + RMS=SQRT(RMS/NUM) + ELSE + RMS=0.0 + ENDIF + + CALL DECLAR('RESULT','DP',' ',DUMMY2,RMS) + CALL DECLAR('VIOLATIONS','DP',' ',DUMMY2,VIOLS) + + RETURN + END diff -u -N -r xplor-nih-2.9.2/source/xangle.fcm PARAxplor-nih-2.9.2/source/xangle.fcm --- xplor-nih-2.9.2/source/xangle.fcm 1970-01-01 01:00:00.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/xangle.fcm 2003-12-05 14:34:36.000000000 +0100 @@ -0,0 +1,59 @@ +C ------------------------------------------------------------------ +C RDC-angles stuff +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + INTEGER MAXXANGLESCLASSES + PARAMETER (MAXXANGLESCLASSES = 40) +C ------------------------------------------------------------------ +C Arrays that hold RDC-angles info +C XANGLESASSNDX tells ending index of the RDC-angles arrays for each +C class +C XANGLESFORCES holds K1 for each class +C ------------------------------------------------------------------ + INTEGER XANGLESASSNDX(MAXXANGLESCLASSES) + REAL XANGLESFORCES(MAXXANGLESCLASSES) + CHARACTER*8 XANGLESCLASSNAMES(MAXXANGLESCLASSES) + LOGICAL PRINTCLASS(MAXXANGLESCLASSES) +C ------------------------------------------------------------------ +C MAXXANGLES is the number of slots set aside for RDC-angles +C assignments +C ANGLESNUM is the total number of RDC-angles entered +C ------------------------------------------------------------------ + INTEGER MAXXANGLES,ANGLESNUM,NCLASSES,CURCLASS,XANGLESMX +C ------------------------------------------------------------------ +C Pointers to arrays to hold atom numbers, observed THETA and PHI, +C and errors +C ------------------------------------------------------------------ + INTEGER XANGLESJPTR,XANGLESKPTR,XANGLESJLST,XANGLESKLST, + & XANGLESOBS1PTR,XANGLESOBS2PTR,XANGLESERRPTR,CALCXANGLESPTR +C ------------------------------------------------------------------ +C Input modes +C ------------------------------------------------------------------ + INTEGER MODE,NEW,UPDATE + PARAMETER (NEW=1) + PARAMETER (UPDATE=2) +C ------------------------------------------------------------------ +C Parameters as set up in ETOR - Not used indeed +C ------------------------------------------------------------------ + DOUBLE PRECISION MCONST + PARAMETER(MCONST=0.0001D0) + DOUBLE PRECISION EPS + PARAMETER(EPS=0.1D0) +C ------------------------------------------------------------------ +C Common blocks +C ------------------------------------------------------------------ + COMMON /CXANGLES/XANGLESCLASSNAMES + COMMON /IXANGLES1/XANGLESASSNDX,MAXXANGLES,ANGLESNUM, + & CURCLASS,NCLASSES,XANGLESMX + COMMON /IXANGLES2/XANGLESJPTR,XANGLESKPTR,XANGLESJLST,XANGLESKLST, + & XANGLESOBS1PTR,XANGLESOBS2PTR, + & XANGLESERRPTR,CALCXANGLESPTR,MODE + COMMON /RXANGLES/XANGLESFORCES + COMMON /LXANGLES/PRINTCLASS + + SAVE /CXANGLES/ + SAVE /IXANGLES1/ + SAVE /IXANGLES2/ + SAVE /RXANGLES/ + SAVE /LXANGLES/ diff -u -N -r xplor-nih-2.9.2/source/xccr.f PARAxplor-nih-2.9.2/source/xccr.f --- xplor-nih-2.9.2/source/xccr.f 1970-01-01 01:00:00.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/xccr.f 2003-12-05 14:34:36.000000000 +0100 @@ -0,0 +1,901 @@ +C ------------------------------------------------------------------ + SUBROUTINE EXCCR (ECCR, WHICH) +C ------------------------------------------------------------------ +C Calls EXCCR2, which does the actual energy calculation +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "xccr.fcm" + INCLUDE "heap.fcm" + + DOUBLE PRECISION ECCR + CHARACTER*7 WHICH + + CALL EXCCR2(ECCR,HEAP(XCCRIPTR),HEAP(XCCRJPTR),HEAP(XCCRKPTR), + & HEAP(XCCRILST),HEAP(XCCRJLST),HEAP(XCCRKLST), + & HEAP(XCCROBSPTR),HEAP(XCCRERRPTR), + & HEAP(CALCXCCRPTR),WHICH) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE EXCCR2 (ECCR,ATOMIPTR,ATOMJPTR,ATOMKPTR, + & ATOMILST,ATOMJLST,ATOMKLST, + & XCCROBS,XCCRERR, + & XCCRCALC,WHICH) +C ------------------------------------------------------------------ +C Calculates cross correlation rate energies +C +C Energies are of the form: +C E = K*(DELTACCR**2) +C where: +C K = FORCE CONSTANT +C DELTACCR = CALCULATED CCR - OBSERVED CCR +C and cross correlation rate function is defined as: +C CCR = KM*(3*COS(THETA)^2-1)/DIST^3 +C +C Note that: +C Atom I = Metal +C Atom J = Detected nucleus (e.g. H) +C Atom K = Coupled nucleus (e.g. N) +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "psf.fcm" + INCLUDE "coord.fcm" + INCLUDE "numbers.fcm" + INCLUDE "deriv.fcm" + INCLUDE "xccr.fcm" + INCLUDE "consta.fcm" + + INTEGER ATOMIPTR(*),ATOMJPTR(*),ATOMKPTR(*) + INTEGER ATOMILST(*),ATOMJLST(*),ATOMKLST(*) + DOUBLE PRECISION XCCROBS(*),XCCRERR(*) + DOUBLE PRECISION XCCRCALC(*),ECCR + CHARACTER*7 WHICH +C ------------------------------------------------------------------ +C Local variables +C ------------------------------------------------------------------ + INTEGER COUNT,CLASS,MM,NN,II,I,J,K,L,NA,M,N + PARAMETER (NA=20) + INTEGER DMM,DNN,OUTFLAG(NA),MCOUNT,NCOUNT,PCOUNT + DOUBLE PRECISION XI,XJ,XK,YI,YJ,YK,ZI,ZJ,ZK,CALCXCCR,E, + & DFXI(NA),DFYI(NA),DFZI(NA),DFXJ(NA),DFYJ(NA), + & DFZJ(NA),DFXK(NA),DFYK(NA),DFZK(NA), + & DJK2,DJI2,DJI5,DJI,DF1,PS,DJK,DJI3, + & O1,O2,O3,O4,O5,O6,O7,O8,O9,O10,O11,O12,O13,O14, + & O15,O16,O17,O18,O19,O20,O21,O22,O23,O24,O25,O26, + & O27,O28,O29,O30,O31,O32,O33,O34,O35,O36,O37,O38, + * O39 + DOUBLE PRECISION OBSXCCR,ERRXCCR,K1,COEF1,A,B,DELTA, + & DP,DELTAXCCR,DD +C ------------------------------------------------------------------ +C Zero out partial energy +C ------------------------------------------------------------------ + ECCR = ZERO + + CLASS = 1 + K1=XCCRFORCES(1) + COEF1=XCCRCOEF1(1) + + COUNT=0 + DO WHILE(COUNT.LT.CCRNUM) + COUNT=COUNT+1 +C ------------------------------------------------------------------ +C Reset individual E to zero +C ------------------------------------------------------------------ + E=0 + + DO WHILE ((XCCRASSNDX(CLASS).LT.COUNT).OR. + & (XCCRASSNDX(CLASS).EQ.0)) + CLASS=CLASS+1 + END DO + + IF (XCCRASSNDX(CLASS).EQ.XCCRASSNDX(CLASS-1)) THEN + COUNT=COUNT-1 + ENDIF + + K1=XCCRFORCES(CLASS) + COEF1=XCCRCOEF1(CLASS) +C ------------------------------------------------------------------ +C Note there should only be one atom for I, J, and K +C ------------------------------------------------------------------ + I=ATOMIPTR(COUNT)+1 + J=ATOMJPTR(COUNT)+1 + K=ATOMKPTR(COUNT)+1 + + XI=X(ATOMILST(I)) + XJ=X(ATOMJLST(J)) + XK=X(ATOMKLST(K)) + YI=Y(ATOMILST(I)) + YJ=Y(ATOMJLST(J)) + YK=Y(ATOMKLST(K)) + ZI=Z(ATOMILST(I)) + ZJ=Z(ATOMJLST(J)) + ZK=Z(ATOMKLST(K)) + + OBSXCCR=XCCROBS(COUNT) + ERRXCCR=XCCRERR(COUNT) +C ------------------------------------------------------------------ +C Initialize calculated cross correlation rates and counter +C ------------------------------------------------------------------ + CALCXCCR=0 + II=0 + MCOUNT=0 + NCOUNT=0 +C ------------------------------------------------------------------ +C Check for correct permutations of paired atoms to get the metal- +C nucleus vector. This depends solely on the order of the assigned +C atoms. OUTFLAG=1 indicates that the permutation is not allowed. +C ------------------------------------------------------------------ + DMM=ATOMJPTR(COUNT+1)-ATOMJPTR(COUNT) + DNN=ATOMKPTR(COUNT+1)-ATOMKPTR(COUNT) + DO MM=ATOMJPTR(COUNT)+1,ATOMJPTR(COUNT+1) + MCOUNT=MCOUNT+1 + NCOUNT=0 + DO NN=ATOMKPTR(COUNT)+1,ATOMKPTR(COUNT+1) + NCOUNT=NCOUNT+1 + II=II+1 + IF ((DMM.GT.1).AND.(DNN.GT.1)) THEN + IF (MCOUNT.EQ.NCOUNT) THEN + OUTFLAG(II)=0 + ELSE + OUTFLAG(II)=1 + ENDIF + + ELSE + OUTFLAG(II)=0 + END IF + END DO + END DO + + II=0 + PCOUNT=0 + DO MM=ATOMJPTR(COUNT)+1,ATOMJPTR(COUNT+1) + DO NN=ATOMKPTR(COUNT)+1,ATOMKPTR(COUNT+1) + II=II+1 + IF (OUTFLAG(II).NE.1) THEN + PCOUNT=PCOUNT+1 +C ------------------------------------------------------------------ +C Calculate the function and its derivatives +C ------------------------------------------------------------------ + O1=XJ-XK + O2=XJ-XI + O3=YJ-YK + O4=YJ-YI + O5=ZJ-ZK + O6=ZJ-ZI + O7=O2**2 + O8=O4**2 + O9=O6**2 + O10=O1**2 + O11=O3**2 + O12=O5**2 + PS=O1*O2+O3*O4+O5*O6 + DJK2=O10+O11+O12 + DJI2=O7+O8+O9 + DJI5=SQRT(DJI2**5) + DJI3=SQRT(DJI2**3) + DJI=SQRT(DJI2) + O13=2.*XJ-XI-XK + O14=6.*PS*O13 + O15=2.*O1*DJI2 + O16=2.*O2*DJK2 + O17=DJK2*DJI5 + O18=3.*PS**2-DJK2*DJI2 + O19=2.*O1*DJI5 + O20=5.*O2*DJK2*DJI3 + O21=O17**2 + O22=6.*PS*(-O2) + O23=6.*PS*(-O1) + O24=2.*YJ-YI-YK + O25=6.*PS*O24 + O26=2.*O3*DJI2 + O27=2.*O4*DJK2 + O28=2.*O3*DJI5 + O29=5.*O4*DJK2*DJI3 + O30=6.*PS*(-O4) + O31=6.*PS*(-O3) + O32=2.*ZJ-ZI-ZK + O33=6.*PS*O32 + O34=2.*O5*DJI2 + O35=2.*O6*DJK2 + O36=2.*O5*DJI5 + O37=5.*O6*DJK2*DJI3 + O38=6.*PS*(-O6) + O39=6.*PS*(-O5) + DJK=SQRT(DJK2) + + DFXI(II)=((O23+O16)*O17-O18*O20)/O21 + DFYI(II)=((O31+O27)*O17-O18*O29)/O21 + DFZI(II)=((O39+O35)*O17-O18*O37)/O21 + DFXJ(II)=((O14-O15-O16)*O17-O18*(O19+O20))/O21 + DFYJ(II)=((O25-O26-O27)*O17-O18*(O28+O29))/O21 + DFZJ(II)=((O33-O34-O35)*O17-O18*(O36+O37))/O21 + DFXK(II)=((O22+O15)*O17+O18*O19)/O21 + DFYK(II)=((O30+O26)*O17+O18*O28)/O21 + DFZK(II)=((O38+O34)*O17+O18*O36)/O21 +C ------------------------------------------------------------------ +C Coefficient +C ------------------------------------------------------------------ + A=COEF1 +C ------------------------------------------------------------------ +C Calculate the cross correlation rate +C ------------------------------------------------------------------ + CALCXCCR=A*((3.*(PS/(DJK*DJI))**2-1.)/DJI3) + END IF + END DO + END DO + + IF (WHICH.EQ.'ANALYZE') THEN + XCCRCALC(COUNT)=CALCXCCR + END IF + + DELTA=(CALCXCCR-OBSXCCR) +C ------------------------------------------------------------------ +C Adjust the deviation based on the error range +C ------------------------------------------------------------------ + IF ((DELTA.LT.0.000).AND.(ABS(DELTA).GT.ERRXCCR)) THEN + DELTAXCCR=DELTA+ERRXCCR + ELSE IF ((DELTA.GT.0.000).AND.(DELTA.GT.ERRXCCR)) THEN + DELTAXCCR=DELTA-ERRXCCR + ELSE + DELTAXCCR=0.0 + END IF + + DD=DELTAXCCR + + II=0 + PCOUNT=0 + DO MM=ATOMJPTR(COUNT)+1,ATOMJPTR(COUNT+1) + DO NN=ATOMKPTR(COUNT)+1,ATOMKPTR(COUNT+1) + II=II+1 + IF (OUTFLAG(II).NE.1) THEN + PCOUNT=PCOUNT+1 + IF (DD.EQ.0.0) THEN + E=E+0.0 + DF1=0.0 + ELSE +C ------------------------------------------------------------------ +C Depending on CCRFLAG, the weight may be proportional to DIST^3 +C ------------------------------------------------------------------ + IF (CCRFLG) THEN + E=E+K1*(DELTAXCCR**2)*DJI3 + DF1=2*K1*DELTAXCCR*DJI3 + ELSE + E=E+K1*(DELTAXCCR**2) + DF1=2*K1*DELTAXCCR + END IF + END IF + END IF + END DO + END DO +C ------------------------------------------------------------------ +C Accumulate energy +C ------------------------------------------------------------------ + ECCR=ECCR+E +C ------------------------------------------------------------------ +C Now update forces if in energy/force mode +C ------------------------------------------------------------------ + IF (WHICH.NE.'ANALYZE') THEN + II=0 + DO MM=ATOMJPTR(COUNT)+1,ATOMJPTR(COUNT+1) + DO NN=ATOMKPTR(COUNT)+1,ATOMKPTR(COUNT+1) + II=II+1 + IF (OUTFLAG(II).NE.1) THEN + DX(ATOMILST(I))=DX(ATOMILST(I))+DF1*DFXI(II) + DY(ATOMILST(I))=DY(ATOMILST(I))+DF1*DFYI(II) + DZ(ATOMILST(I))=DZ(ATOMILST(I))+DF1*DFZI(II) + DX(ATOMJLST(J))=DX(ATOMJLST(J))+DF1*DFXJ(II) + DY(ATOMJLST(J))=DY(ATOMJLST(J))+DF1*DFYJ(II) + DZ(ATOMJLST(J))=DZ(ATOMJLST(J))+DF1*DFZJ(II) + DX(ATOMKLST(K))=DX(ATOMKLST(K))+DF1*DFXK(II) + DY(ATOMKLST(K))=DY(ATOMKLST(K))+DF1*DFYK(II) + DZ(ATOMKLST(K))=DZ(ATOMKLST(K))+DF1*DFZK(II) + END IF + END DO + END DO + END IF + END DO + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE READXCCR +C ------------------------------------------------------------------ +C Reads in cross correlation rate information +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "comand.fcm" + INCLUDE "xccr.fcm" + INCLUDE "funct.fcm" + INCLUDE "psf.fcm" + INCLUDE "heap.fcm" + INCLUDE "numbers.fcm" + + INTEGER COUNT,SPTR,OLDCLASS,OLDMAXXCCR,FNCLASS + DOUBLE PRECISION K1,CUTOFF,COEF1 + CHARACTER*4 THENAME + CHARACTER*132 NOMEFILE + INTEGER WEISWI + + NOMEFILE='XCCR.OUT' + + SPTR=ALLHP(INTEG4(NATOM)) + CALL PUSEND('XCCR>') +862 CONTINUE + CALL NEXTWD('XCCR>') + CALL MISCOM('XCCR>',USED) + IF (.NOT.USED) THEN +C ------------------------------------------------------------------ +C Documentation +C ------------------------------------------------------------------ + IF (WD(1:4).EQ.'HELP') THEN + WRITE(DUNIT,'(10X,A)') + & ' XCCR {} END ', + & ' :== ', + & ' ASSIGN ', + & ' * Restraint: Metal Coupled_nucleus', + & ' Detected_nucleus CCR Err *', + & ' CLASsification ', + & ' * Starts a new class *', + & ' WEIP <1|0>', + & ' * Switch weight proportional to R^3 (1=Y,0=N) *', + & ' COEFficient ', + & ' * Proportionality constant *', + & ' FORCeconstant ', + & ' * Force constant for the current class *', + & ' NREStraints ', + & ' * Number of slots for restraints to allocate *', + & ' PRINt THREshold ', + & ' * Prints violations larger than the THREshold *', + & ' RESEt', + & ' * Erases the restraint table, but keeps NRES *', + & ' FRUN', + & ' * Runs FANTACROSS *' +C ------------------------------------------------------------------ +C About FANTACROSS +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'FRUN') THEN + CALL NEXTI('FANTACROSS on class number =',FNCLASS) + IF (FNCLASS.GT.NCLASSES.OR.FNCLASS.EQ.0) THEN + PRINT*,'%FRUN-ERR: This class does not exist...' + ELSE + CALL FANTACCR(HEAP(XCCRIPTR),HEAP(XCCRJPTR), + & HEAP(XCCRKPTR),HEAP(XCCRILST), + & HEAP(XCCRJLST),HEAP(XCCRKLST), + & HEAP(XCCROBSPTR),HEAP(XCCRERRPTR), + & FNCLASS,NOMEFILE) + END IF +C ------------------------------------------------------------------ +C Get class name. Determine if it's an already-defined class. +C Insert a new class if it's not. +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'CLAS') THEN + OLDCLASS=CURCLASS + CALL NEXTA4('Class name =',THENAME) + MODE=NEW + DO COUNT=1,NCLASSES + IF (XCCRCLASSNAMES(COUNT).EQ.THENAME) THEN + MODE=UPDATE + CURCLASS=COUNT + END IF + END DO + IF (MODE.EQ.NEW) THEN +C ------------------------------------------------------------------ +C Make sure you can't add more than the maximum number of classes +C ------------------------------------------------------------------ + IF (OLDCLASS.EQ.MAXXCCRCLASSES) THEN + CALL DSPERR('XCCR','Too many classes.') + CALL DSPERR('XCCR', + & 'Increase MAXXCCRCLASSES and recompile.') + CALL WRNDIE(-5, 'READXCCR', + & 'Too many CCR classes.') + END IF + NCLASSES=NCLASSES+1 + CURCLASS=NCLASSES + XCCRCLASSNAMES(CURCLASS)=THENAME +C ------------------------------------------------------------------ +C If this isn't the first class, close off the old class +C ------------------------------------------------------------------ + IF (NCLASSES.GT.1) THEN + XCCRASSNDX(OLDCLASS)=CCRNUM + END IF + END IF +C ------------------------------------------------------------------ +C Set force constant for current class +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'FORC') THEN +C ------------------------------------------------------------------ +C Start a default class if there isn't one defined +C ------------------------------------------------------------------ + IF (CURCLASS.EQ.0) THEN + NCLASSES = 1 + CURCLASS = 1 + END IF + CALL NEXTF('Force constant =',K1) + WRITE(DUNIT,'(A,A,A,F8.3)') + & 'Setting force const for class ', + & XCCRCLASSNAMES(CURCLASS),' to ',K1 + XCCRFORCES(CURCLASS)=K1 +C ------------------------------------------------------------------ +C Set weight proportional to DIST^3 +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'WEIP') THEN + CALL NEXTI('Weight proportional to R^3 (Y=1,N=0):',WEISWI) + IF (WEISWI.EQ.1) THEN + CCRFLG=.TRUE. + WRITE (DUNIT,'(A)') 'CCR weight proportional to R^3.' + ELSE IF (WEISWI.EQ.0) THEN + CCRFLG=.FALSE. + WRITE (DUNIT,'(A)') 'CCR weight not dependent on R.' + ELSE + WRITE (DUNIT,'(A)') 'Unknown switch. Default weight.' + END IF +C ------------------------------------------------------------------ +C Set coefficient constant for current class +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'COEF') THEN + CALL NEXTF('CCR COEFficient =',COEF1) +C ------------------------------------------------------------------ +C Start a default class if there isn't one already defined +C ------------------------------------------------------------------ + IF (CURCLASS.EQ.0) THEN + NCLASSES = 1 + CURCLASS = 1 + END IF + WRITE(DUNIT,'(A,A,A,F8.3)') + & 'Setting coefficient for class ', + & XCCRCLASSNAMES(CURCLASS),'to',COEF1 + XCCRCOEF1(CURCLASS)=COEF1 +C ------------------------------------------------------------------ +C Reset cross correlation rates database +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'RESE') THEN + CALL XCCRDEFAULTS + CALL ALLOCXCCR(0,MAXXCCR) +C ------------------------------------------------------------------ +C Change number of assignment slots +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'NRES') THEN + OLDMAXXCCR=MAXXCCR + CALL NEXTI('Number of slots =',MAXXCCR) + CALL ALLOCXCCR(OLDMAXXCCR,MAXXCCR) + WRITE(DUNIT,'(A,I8,A)') + & 'XCCR: Allocating space for',MAXXCCR, + & 'number of restraints.' +C ------------------------------------------------------------------ +C Read in an assignment +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'ASSI') THEN +C ------------------------------------------------------------------ +C Make sure you can't add more assignments than you have slots for +C ------------------------------------------------------------------ + IF (XCCRMX.EQ.MAXXCCR) THEN + CALL DSPERR('XCCR','Too many assignments.') + CALL DSPERR('XCCR', + & 'Increasing NREStraints by 100.') + OLDMAXXCCR=MAXXCCR + MAXXCCR=MAXXCCR+100 + CALL ALLOCXCCR(OLDMAXXCCR,MAXXCCR) + END IF +C ------------------------------------------------------------------ +C If there isn't a class specified, start a default class +C ------------------------------------------------------------------ + IF (CURCLASS.EQ.0) THEN + NCLASSES=1 + CURCLASS=1 + END IF + CALL READXCCR2(HEAP(XCCRIPTR),HEAP(XCCRJPTR), + & HEAP(XCCRKPTR),HEAP(XCCRILST), + & HEAP(XCCRJLST),HEAP(XCCRKLST), + & HEAP(XCCROBSPTR),HEAP(XCCRERRPTR), + & HEAP(SPTR)) +C ------------------------------------------------------------------ +C Print violations +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'PRIN') THEN + CALL NEXTWD('PRINT>') + IF (WD(1:4).NE.'THRE') THEN + CALL DSPERR('XCCR', + & 'PRINt expects THREshold parameter.') + ELSE + CALL NEXTF('THREshold =',CUTOFF) + IF (CUTOFF.LT.ZERO) THEN + CALL DSPERR('XCCR', + & 'Cutoff must be positive.') + CUTOFF = ABS(CUTOFF) + END IF + CALL NEXTA4('ALL OR CLASs>',THENAME) + IF (THENAME(1:3).EQ.'ALL') THEN + DO COUNT=1,NCLASSES + PRINTCLASS(COUNT)=.TRUE. + END DO + ELSE IF (THENAME(1:4).EQ.'CLAS') THEN + CALL NEXTA4('Class name =',THENAME) + DO COUNT=1,NCLASSES + IF (XCCRCLASSNAMES(COUNT).EQ.THENAME) THEN + PRINTCLASS(COUNT)=.TRUE. + ELSE + PRINTCLASS(COUNT)=.FALSE. + END IF + END DO + ELSE + DO COUNT=1,NCLASSES + PRINTCLASS(COUNT)=.TRUE. + END DO + END IF + CALL PRINTXCCR(CUTOFF,HEAP(CALCXCCRPTR), + & HEAP(XCCROBSPTR),HEAP(XCCRERRPTR), + & HEAP(XCCRIPTR),HEAP(XCCRJPTR), + & HEAP(XCCRKPTR),HEAP(XCCRILST), + & HEAP(XCCRJLST),HEAP(XCCRKLST)) + END IF +C ------------------------------------------------------------------ +C Check for END statement +C ------------------------------------------------------------------ + ELSE + CALL CHKEND('XCCR>',DONE) + END IF + END IF + IF (.NOT.(DONE)) GOTO 862 + DONE=.FALSE. + CALL FREHP(SPTR,INTEG4(NATOM)) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE ALLOCXCCR (OLDSIZE,NEWSIZE) +C ------------------------------------------------------------------ +C Resets cross correlation rate arrays to hold size entries +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "funct.fcm" + INCLUDE "xccr.fcm" + + INTEGER OLDSIZE,NEWSIZE + + IF (OLDSIZE.NE.0) THEN + CALL FREHP(XCCRIPTR,INTEG4(OLDSIZE)) + CALL FREHP(XCCRJPTR,INTEG4(OLDSIZE)) + CALL FREHP(XCCRKPTR,INTEG4(OLDSIZE)) + + CALL FREHP(XCCRILST,INTEG4(OLDSIZE)) + CALL FREHP(XCCRJLST,INTEG4(OLDSIZE)) + CALL FREHP(XCCRKLST,INTEG4(OLDSIZE)) + + CALL FREHP(XCCROBSPTR,IREAL8(OLDSIZE)) + CALL FREHP(XCCRERRPTR,IREAL8(OLDSIZE)) + CALL FREHP(CALCXCCRPTR,IREAL8(OLDSIZE)) + END IF + XCCRIPTR=ALLHP(INTEG4(NEWSIZE)) + XCCRJPTR=ALLHP(INTEG4(NEWSIZE)) + XCCRKPTR=ALLHP(INTEG4(NEWSIZE)) + + XCCRILST=ALLHP(INTEG4(NEWSIZE)) + XCCRJLST=ALLHP(INTEG4(NEWSIZE)) + XCCRKLST=ALLHP(INTEG4(NEWSIZE)) + + XCCROBSPTR=ALLHP(IREAL8(NEWSIZE)) + XCCRERRPTR=ALLHP(IREAL8(NEWSIZE)) + CALCXCCRPTR=ALLHP(IREAL8(NEWSIZE)) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE XCCRDEFAULTS +C ------------------------------------------------------------------ +C Sets up defaults +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "xccr.fcm" + + INTEGER COUNT + + MODE=NEW + MAXXCCR=200 + XCCRMX=200 + CCRNUM=0 + NCLASSES=0 + CURCLASS=0 + CCRFLG=.FALSE. + DO COUNT=1,MAXXCCRCLASSES + XCCRCLASSNAMES(COUNT)='DEFAULT' + XCCRASSNDX(COUNT)=0 + XCCRFORCES(COUNT)=5.0 + XCCRCOEF1(COUNT)=1.0 + END DO + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE READXCCR2(ATOMIPTR,ATOMJPTR,ATOMKPTR, + & ATOMILST,ATOMJLST,ATOMKLST, + & XCCROBS,XCCRERR,SEL) +C ------------------------------------------------------------------ +C Reads actual cross correlation rate assignments into arrays +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "psf.fcm" + INCLUDE "coord.fcm" + INCLUDE "xccr.fcm" + + INTEGER ATOMIPTR(*),ATOMJPTR(*),ATOMKPTR(*) + INTEGER ATOMILST(*),ATOMJLST(*),ATOMKLST(*) + INTEGER SEL(*) + INTEGER ITMP,JTMP,KTMP,II + DOUBLE PRECISION XCCROBS(*),XCCRERR(*) +C ------------------------------------------------------------------ +C Local variables +C ------------------------------------------------------------------ + INTEGER NSEL,INSERTPOS,COUNT,CURSTOP,OTHERSTOP,NFLAG + INTEGER I + DOUBLE PRECISION XCCRO,XCCRE +C ------------------------------------------------------------------ +C If we're in UPDATE mode, make a space for the new line +C ------------------------------------------------------------------ + NFLAG = 0 + IF (MODE.EQ.UPDATE) THEN + DO COUNT=CCRNUM+1,XCCRASSNDX(CURCLASS)+1,-1 + ATOMIPTR(COUNT)=ATOMIPTR(COUNT-1) + ATOMJPTR(COUNT)=ATOMJPTR(COUNT-1) + ATOMKPTR(COUNT)=ATOMKPTR(COUNT-1) + XCCROBS(COUNT)=XCCROBS(COUNT-1) + XCCRERR(COUNT)=XCCRERR(COUNT-1) + END DO + CURSTOP=XCCRASSNDX(CURCLASS) + DO COUNT=1,NCLASSES + OTHERSTOP=XCCRASSNDX(COUNT) + IF (OTHERSTOP.GT.CURSTOP) THEN + XCCRASSNDX(COUNT)=OTHERSTOP+1 + END IF + END DO + XCCRASSNDX(CURCLASS)=CURSTOP+1 + INSERTPOS=CURSTOP + CCRNUM=CCRNUM+1 + ELSE + CCRNUM=CCRNUM+1 + INSERTPOS=CCRNUM + XCCRASSNDX(CURCLASS)=INSERTPOS + END IF +C ------------------------------------------------------------------ +C Reading in the atom selection in the restraint table +C ------------------------------------------------------------------ + CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.) + IF (NSEL.GT.1) THEN + CALL DSPERR('XCCR', + & 'More than 1 atom in SEL for atom I. Using first') + END IF + IF (NSEL.EQ.0) THEN + CALL DSPERR('XCCR','Error with atom selection') + NFLAG = 1 + END IF + CALL MAKIND(SEL, NATOM, NSEL) + IF (INSERTPOS.EQ.1) ATOMIPTR(INSERTPOS)=0 + II=ATOMIPTR(INSERTPOS) + II=II+1 + IF (II.GT.XCCRMX) XCCRMX=II + ATOMILST(II) = SEL(1) + ATOMIPTR(INSERTPOS+1)=II + + CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.) + IF (NSEL.GT.1) THEN + CALL DSPERR('XCCR', + & 'More than 1 atom in SEL for atom J. Using first') + END IF + IF (NSEL.EQ.0) THEN + CALL DSPERR('XCCR','Error with atom selection') + NFLAG = 1 + END IF + CALL MAKIND(SEL,NATOM,NSEL) + IF (INSERTPOS.EQ.1) ATOMJPTR(INSERTPOS)=0 + II=ATOMJPTR(INSERTPOS) + II=II+1 + IF (II.GT.XCCRMX) XCCRMX=II + ATOMJLST(II)=SEL(1) + ATOMJPTR(INSERTPOS+1)=II + + CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.) + IF (NSEL.GT.1) THEN + CALL DSPERR('XCCR', + & 'More than 1 atom in SEL for atom K. Using first') + END IF + IF (NSEL.EQ.0) THEN + CALL DSPERR('XCCR','Error with atom selection') + NFLAG=1 + END IF + CALL MAKIND(SEL,NATOM,NSEL) + IF (INSERTPOS.EQ.1) ATOMKPTR(INSERTPOS)=0 + II=ATOMKPTR(INSERTPOS) + II=II+1 + IF (II.GT.XCCRMX) XCCRMX=II + ATOMKLST(II)=SEL(1) + ATOMKPTR(INSERTPOS+1)=II +C ------------------------------------------------------------------ +C Reading in the observed cross correlation rate +C ------------------------------------------------------------------ + CALL NEXTF('Observed CCR =',XCCRO) + CALL NEXTF('Error in CCR =',XCCRE) + + XCCROBS(INSERTPOS)=XCCRO + XCCRERR(INSERTPOS)=XCCRE +C ------------------------------------------------------------------ +C Check for error in atom selection. If there is one then reset the +C counter for restraint +C ------------------------------------------------------------------ + IF (NFLAG.EQ.1) THEN + CCRNUM=CCRNUM-1 + END IF + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE XCCRINIT +C ------------------------------------------------------------------ +C Initializes cross correlation rate stuff +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "xccr.fcm" + + CALL XCCRDEFAULTS + CALL ALLOCXCCR(0,MAXXCCR) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE PRINTXCCR (CUTOFF,XCCRCALC,XCCROBS,XCCRERR, + & ATOMIPTR,ATOMJPTR,ATOMKPTR, + & ATOMILST,ATOMJLST,ATOMKLST) +C ------------------------------------------------------------------ +C Prints cross correlation rates with dveiation larger than cutoff, +C calculates RMSD and puts it into $RESULT +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "xccr.fcm" + INCLUDE "comand.fcm" + INCLUDE "psf.fcm" + INCLUDE "numbers.fcm" + + DOUBLE PRECISION CUTOFF,XCCRCALC(*),XCCROBS(*),XCCRERR(*) + INTEGER ATOMILST(*),ATOMJLST(*),ATOMKLST(*) + INTEGER ATOMIPTR(*),ATOMJPTR(*),ATOMKPTR(*) +C ------------------------------------------------------------------ +C Local variables +C ------------------------------------------------------------------ + DOUBLE PRECISION CALCXCCR,OBSXCCR,DELTAXCCR,DELTA,DP + INTEGER COUNT,CLASS,I,J,K,NUM,II + DOUBLE PRECISION RMS,VIOLS,ERRXCCR,XCCRENERGY + DOUBLE COMPLEX DUMMY2 + LOGICAL PRINTTHISCLASS + + RMS=ZERO + VIOLS=ZERO + NUM=0 +C ------------------------------------------------------------------ +C Make sure that the array of calculated CCR is up to date +C ------------------------------------------------------------------ + CALL EXCCR(XCCRENERGY,'ANALYZE') + WRITE (PUNIT,'(A)') 'The following cross correlation rates have' + WRITE (PUNIT,'(A)') 'deviation larger than or' + WRITE (PUNIT,'(A)') 'equal to the cutoff:' +C ------------------------------------------------------------------ +C Write out first class heading +C ------------------------------------------------------------------ + CLASS=1 + PRINTTHISCLASS=PRINTCLASS(CLASS) + IF (PRINTTHISCLASS) THEN + WRITE (PUNIT,'(A,A)') 'Class ',XCCRCLASSNAMES(1) + END IF +C ------------------------------------------------------------------ +C For every cross correlation rate entry... +C ------------------------------------------------------------------ + COUNT = 0 + DO WHILE(COUNT.LT.CCRNUM) + COUNT=COUNT+1 +C ------------------------------------------------------------------ +C Is this the start of a new class? +C ------------------------------------------------------------------ + IF (XCCRASSNDX(CLASS).LT.COUNT) THEN + CLASS=CLASS+1 + PRINTTHISCLASS=PRINTCLASS(CLASS) + IF (XCCRASSNDX(CLASS).EQ.XCCRASSNDX(CLASS-1)) THEN + COUNT=COUNT - 1 + END IF + IF (PRINTTHISCLASS) THEN + WRITE (PUNIT,'(A,A)') 'Class ',XCCRCLASSNAMES(CLASS) + END IF + END IF +C ------------------------------------------------------------------ +C If this assignment is in a class to be printed +C and make sure there is an entry for that class +C ------------------------------------------------------------------ + IF ((PRINTTHISCLASS).AND. + & (XCCRASSNDX(CLASS).NE.XCCRASSNDX(CLASS-1))) THEN +C ------------------------------------------------------------------ +C Always update RMSD +C ------------------------------------------------------------------ + CALCXCCR=XCCRCALC(COUNT) + OBSXCCR=XCCROBS(COUNT) + ERRXCCR=XCCRERR(COUNT) + DP=(CALCXCCR-OBSXCCR) + + IF ((DP.LT.0.000).AND.(ABS(DP).GT.ERRXCCR)) THEN + DELTAXCCR=DP+ERRXCCR + ELSE IF ((DP.GT.0.000).AND.(DP.GT.ERRXCCR)) THEN + DELTAXCCR=DP-ERRXCCR + ELSE + DELTAXCCR=0.0 + END IF + + RMS=RMS+DELTAXCCR**2 + NUM=NUM+1 +C ------------------------------------------------------------------ +C Print out deviations larger than cutoff +C and update number of violations +C ------------------------------------------------------------------ + IF (ABS(DELTAXCCR).GE.CUTOFF) THEN + I=ATOMILST(ATOMIPTR(COUNT)+1) + J=ATOMJLST(ATOMJPTR(COUNT)+1) + K=ATOMKLST(ATOMKPTR(COUNT)+1) + WRITE(PUNIT,'(A,A)') '==============================', + & '==============================' + WRITE(PUNIT,'(A)') ' Set-I-atoms' + DO II=ATOMIPTR(COUNT)+1,ATOMIPTR(COUNT+1) + I=ATOMILST(II) + WRITE(PUNIT,'(9X,4(1X,A))') SEGID(I),RESID(I), + & RES(I),TYPE(I) + END DO + WRITE(PUNIT,'(A)') ' Set-J-atoms' + DO II=ATOMJPTR(COUNT)+1,ATOMJPTR(COUNT+1) + J=ATOMJLST(II) + WRITE(PUNIT,'(9X,4(1X,A))') SEGID(J),RESID(J), + & RES(J),TYPE(J) + END DO + WRITE(PUNIT,'(A)') ' Set-K-atoms' + DO II=ATOMKPTR(COUNT)+1,ATOMKPTR(COUNT+1) + K=ATOMKLST(II) + WRITE(PUNIT,'(9X,4(1X,A))') SEGID(K),RESID(K), + & RES(K),TYPE(K) + END DO + WRITE(PUNIT,'(2(2X,A,1X,F8.3))') + & 'Calc: ',CALCXCCR,'Obs: ',OBSXCCR + WRITE(PUNIT,'(2X,A,1X,F8.3,2X,A,1X,F8.3)') + & 'Error: ',ERRXCCR,'Delta: ',DELTAXCCR + VIOLS=VIOLS+ONE + END IF + END IF + END DO + + IF (NUM.GT.0) THEN + RMS=SQRT(RMS/NUM) + ELSE + RMS=0.0 + ENDIF + + CALL DECLAR('RESULT','DP',' ',DUMMY2,RMS) + CALL DECLAR('VIOLATIONS','DP',' ',DUMMY2,VIOLS) + + RETURN + END diff -u -N -r xplor-nih-2.9.2/source/xccr.fcm PARAxplor-nih-2.9.2/source/xccr.fcm --- xplor-nih-2.9.2/source/xccr.fcm 1970-01-01 01:00:00.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/xccr.fcm 2003-12-05 14:34:36.000000000 +0100 @@ -0,0 +1,58 @@ +C ------------------------------------------------------------------ +C Cross correlation rate stuff +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + INTEGER MAXXCCRCLASSES + PARAMETER (MAXXCCRCLASSES = 40) +C ------------------------------------------------------------------ +C Arrays that hold cross correlation rate info +C XCCRASSNDX tells ending index of the CCR arrays for each class +C XCCRFORCES holds K1 for each class +C ------------------------------------------------------------------ + INTEGER XCCRASSNDX (MAXXCCRCLASSES) + REAL XCCRFORCES(MAXXCCRCLASSES) + REAL XCCRCOEF1(MAXXCCRCLASSES) + CHARACTER*8 XCCRCLASSNAMES(MAXXCCRCLASSES) + LOGICAL PRINTCLASS(MAXXCCRCLASSES),CCRFLG +C ------------------------------------------------------------------ +C MAXXCCR is the number of slots set aside for CCR assignments +C CCRNUM is the total number of CCR entered +C ------------------------------------------------------------------ + INTEGER MAXXCCR,CCRNUM,NCLASSES,CURCLASS,XCCRMX +C ------------------------------------------------------------------ +C Pointers to arrays to hold atom numbers, observed CCR and errors +C ------------------------------------------------------------------ + INTEGER XCCRIPTR,XCCRJPTR,XCCRKPTR,XCCRILST,XCCRJLST,XCCRKLST, + & XCCROBSPTR,XCCRERRPTR,CALCXCCRPTR +C ------------------------------------------------------------------ +C Input modes +C ------------------------------------------------------------------ + INTEGER MODE,NEW,UPDATE + PARAMETER (NEW=1) + PARAMETER (UPDATE=2) +C ------------------------------------------------------------------ +C Parameters as set up in ETOR - Not used indeed +C ------------------------------------------------------------------ + DOUBLE PRECISION MCONST + PARAMETER(MCONST=0.0001D0) + DOUBLE PRECISION EPS + PARAMETER(EPS=0.1D0) +C ------------------------------------------------------------------ +C Common blocks +C ------------------------------------------------------------------ + COMMON /CXCCR/XCCRCLASSNAMES + COMMON /IXCCR1/XCCRASSNDX,MAXXCCR,CCRNUM,CURCLASS,NCLASSES,XCCRMX + COMMON /IXCCR2/XCCRIPTR,XCCRJPTR,XCCRKPTR, + & XCCRILST,XCCRJLST,XCCRKLST, + & XCCROBSPTR,XCCRERRPTR,CALCXCCRPTR,MODE + COMMON /RXCCR/XCCRFORCES,XCCRCOEF1 + COMMON /LXCCR/PRINTCLASS + COMMON /WXCCR/CCRFLG + + SAVE /CXCCR/ + SAVE /IXCCR1/ + SAVE /IXCCR2/ + SAVE /RXCCR/ + SAVE /LXCCR/ + SAVE /WXCCR/ diff -u -N -r xplor-nih-2.9.2/source/xdipo_pcs.f PARAxplor-nih-2.9.2/source/xdipo_pcs.f --- xplor-nih-2.9.2/source/xdipo_pcs.f 1970-01-01 01:00:00.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/xdipo_pcs.f 2003-12-05 14:34:36.000000000 +0100 @@ -0,0 +1,1538 @@ +C ------------------------------------------------------------------ + SUBROUTINE FANTALIN (ATOMIPTR,ATOMJPTR,ATOMKPTR,ATOMLPTR,ATOMMPTR, + & ATOMNPTR,ATOMILST,ATOMJLST,ATOMKLST,ATOMLLST, + & ATOMMLST,ATOMNLST,XPCSOBS,XPCSERR,NCLASS, + & NOMEFILE) +C ------------------------------------------------------------------ +C Fits pseudocontact shift tensor on a given structure +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "psf.fcm" + INCLUDE "coord.fcm" + INCLUDE "numbers.fcm" + INCLUDE "deriv.fcm" + INCLUDE "xdipo_pcs.fcm" + INCLUDE "consta.fcm" + INCLUDE "fantaxplor.fcm" + INCLUDE 'heap.fcm' + + DOUBLE PRECISION APARA,APERP + COMMON /FRUN/APERP,APARA + LOGICAL FSWITCH,FSAVE,FERR + COMMON /FPARAM/FSWITCH,FSAVE,FERR + INTEGER FVIOLS + COMMON /VIOL/FVIOLS + DOUBLE PRECISION FENERGTOT + COMMON /FENERG/ FENERGTOT + DOUBLE PRECISION XPCSE + COMMON /XPCSENERGY/ XPCSE + DOUBLE PRECISION FENT(20,2000),FPER(20,2000),FPAR(20,2000), + & FXPCSE(20,2000) + INTEGER CONTSAVE(20) + COMMON /FMEDIA/FENT,FPER,FPAR,FXPCSE,CONTSAVE + CHARACTER NOMEFILE*132 + INTEGER ATOMIPTR(*),ATOMJPTR(*),ATOMKPTR(*),ATOMLPTR(*) + INTEGER ATOMMPTR(*),ATOMNPTR(*) + INTEGER ATOMILST(*),ATOMJLST(*),ATOMKLST(*),ATOMLLST(*) + INTEGER ATOMMLST(*),ATOMNLST(*) + INTEGER COUNT,MM,NN,I,J,K,L,NA,M,N,CLASS,NORES + INTEGER NCLASS,NNCO,NNCC + DOUBLE PRECISION XPCSOBS(*),XPCSERR(*) + DOUBLE PRECISION OBSXPCS,ERRXPCS + DOUBLE PRECISION XI,XJ,XK,XL,XM,XN,YI,YJ,YK,YL,YM,YN, + & ZI,ZJ,ZK,ZL,ZM,ZN + DOUBLE PRECISION XFA(1000),YFA(1000),ZFA(1000), + & XMFA(1000),YMFA(1000),ZMFA(1000), + & ERRFA(1000),OBSFA(1000), + & WPROT(1000) + DOUBLE COMPLEX DUMMY2 +C ------------------------------------------------------------------ +C CLASS must remain this during the cycle of variables exchange. +C The class is defined by NCLASS. +C ------------------------------------------------------------------ + CLASS=1 + NNCO=0 + NNCC=0 + COUNT=0 + + PRINT*,'This is FANTALIN for PCS.' + PRINT*,'PCSNUM is:',PCSNUM + + DO WHILE (COUNT.LT.PCSNUM) + COUNT=COUNT+1 + DO WHILE ((XPCSASSNDX(CLASS).LT.COUNT).OR. + & (XPCSASSNDX(CLASS).EQ.0)) + CLASS=CLASS+1 + END DO + IF (XPCSASSNDX(CLASS).EQ.XPCSASSNDX(CLASS-1)) THEN + COUNT=COUNT-1 + ENDIF + + IF (NCLASS.EQ.CLASS) THEN + NNCO=NNCO+1 + I=ATOMIPTR(COUNT)+1 + J=ATOMJPTR(COUNT)+1 + K=ATOMKPTR(COUNT)+1 + L=ATOMLPTR(COUNT)+1 + XI=X(ATOMILST(I)) + XJ=X(ATOMJLST(J)) + XK=X(ATOMKLST(K)) + XL=X(ATOMLLST(L)) + YI=Y(ATOMILST(I)) + YJ=Y(ATOMJLST(J)) + YK=Y(ATOMKLST(K)) + YL=Y(ATOMLLST(L)) + ZI=Z(ATOMILST(I)) + ZJ=Z(ATOMJLST(J)) + ZK=Z(ATOMKLST(K)) + ZL=Z(ATOMLLST(L)) + OBSFA(NNCO)=XPCSOBS(COUNT) + OBSXPCS=XPCSOBS(COUNT) + ERRFA(NNCO)=XPCSERR(COUNT) + ERRXPCS=XPCSERR(COUNT) + END IF + + DO MM=ATOMMPTR(COUNT)+1,ATOMMPTR(COUNT+1) + DO NN=ATOMNPTR(COUNT)+1,ATOMNPTR(COUNT+1) + IF (NCLASS.EQ.CLASS) THEN + NNCC=NNCC+1 + XMFA(NNCC)=X(ATOMMLST(MM)) + XFA(NNCC)=X(ATOMNLST(NN)) + YMFA(NNCC)=Y(ATOMMLST(MM)) + YFA(NNCC)=Y(ATOMNLST(NN)) + ZMFA(NNCC)=Z(ATOMMLST(MM)) + ZFA(NNCC)=Z(ATOMNLST(NN)) + NORES=ATOMNLST(NN) + NRESIDFA(NNCC)=RESID(NORES) + NRESFA(NNCC)=RES(NORES) + NTYPEFA(NNCC)=TYPE(NORES) + END IF + END DO + END DO + END DO + + IF (NNCC.EQ.NNCO) THEN + CALL FANTALINEXE(XFA,YFA,ZFA,XMFA,YMFA,ZMFA,NNCO,OBSFA,ERRFA, + & 1,NOMEFILE,WPROT,0,PCSFLG) + PRINT*,'FERR is:',FERR + IF (FERR.EQ.1) THEN + CALL FANTAERRORE(XFA,YFA,ZFA,XMFA,YMFA,ZMFA,NNCO,OBSFA, + & ERRFA,1,PCSFLG) + END IF + ELSE + PRINT*,'Error in atom count...' + END IF + + PRINT *,'FSWITCH is:',FSWITCH + IF (FSWITCH.EQ.1) THEN + XPCSCOEF1(NCLASS)=APARA + XPCSCOEF2(NCLASS)=APERP + PRINT*,'For class:',NCLASS + PRINT*,'New value for A1 =',APARA + PRINT*,'New value for A2 =',APERP + END IF + + IF (FSAVE.EQ.1) THEN + CONTSAVE(NCLASS)=CONTSAVE(NCLASS)+1 + FENT(NCLASS,CONTSAVE(NCLASS))=FENERGTOT + FXPCSE(NCLASS,CONTSAVE(NCLASS))=XPCSE + FPAR(NCLASS,CONTSAVE(NCLASS))=APARA + FPER(NCLASS,CONTSAVE(NCLASS))=APERP + PRINT*,'Value of A1 =',FPAR(NCLASS,CONTSAVE(NCLASS)) + PRINT*,'Value of A2 =',FPER(NCLASS,CONTSAVE(NCLASS)) + PRINT*,'Total energy =',FENERGTOT + PRINT*,'XDIPO_PCS energy =',XPCSE + PRINT*,'CONTSAVE is:',CONTSAVE(NCLASS) + END IF + + CALL DECLAR('CHIAX','DP',' ',DUMMY2,APARA) + CALL DECLAR('CHIRH','DP',' ',DUMMY2,APERP) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE FANTAMEDIA (FMA,PERC) +C ------------------------------------------------------------------ +C Averages tensor parameters on a given sub-ensemble of structures +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + DOUBLE PRECISION VETE(2,2000),VETT(3,2000) + DOUBLE PRECISION FENT(20,2000),FPER(20,2000),FPAR(20,2000), + & FXPCSE(20,2000) + COMMON /FMEDIA/FENT,FPER,FPAR,FXPCSE,CONTSAVE + INTEGER CONTSAVE(20) + INTEGER FMA,COUNT,CB,COUNT1,NUM,CONFPRE + DOUBLE PRECISION PERC + DOUBLE PRECISION M1PAR,M1PER,DEVPAR + DOUBLE COMPLEX DUMMY2 +C ------------------------------------------------------------------ +C Through FMA the class on which the average must be done is passed. +C If it is zero, all the counters are reset. +C ------------------------------------------------------------------ + IF (FMA.EQ.0) THEN + DO COUNT=1,19 + CONTSAVE(COUNT)=0 + END DO + PRINT*,'Mean vectors reset.' + ELSE +C ------------------------------------------------------------------ +C The best are selected +C ------------------------------------------------------------------ + DO COUNT=1,CONTSAVE(FMA) + VETE(1,COUNT)=0 + END DO + + DO COUNT=1,CONTSAVE(FMA) + NUM=CONTSAVE(FMA) + DO WHILE (VETE(1,NUM).EQ.1) + NUM=NUM-1 + END DO + DO COUNT1=1,CONTSAVE(FMA)-1 + IF ((FENT(FMA,COUNT1).LE.FENT(FMA,NUM)).AND. + & (VETE(1,COUNT1).EQ.0.0)) THEN + NUM=COUNT1 + END IF + END DO + PRINT*,'NUM is:',NUM + VETE(1,NUM)=1 + VETE(2,COUNT)=FENT(FMA,NUM) + VETT(1,COUNT)=FPAR(FMA,NUM) + VETT(2,COUNT)=FPER(FMA,NUM) + END DO + + CONFPRE=NINT((DBLE(CONTSAVE(FMA))*(PERC/100.))) + IF (CONFPRE.EQ.0) THEN + CONFPRE=1 + END IF + PRINT*,'CONFPRE is:',CONFPRE + PRINT*,'PERC is:',PERC + PRINT*,'DBLE(CONTSAVE(FMA)) is:',DBLE(CONTSAVE(FMA)) + + DO COUNT=1,CONTSAVE(FMA) + PRINT*,'NUMBER',COUNT,'ENERGY',VETE(1,COUNT),VETE(2,COUNT) + PRINT*,'A1',VETT(1,COUNT),' A2',VETT(2,COUNT) + END DO + + CB=0 + M1PAR=0 + M1PER=0 + DO COUNT=1,CONFPRE + M1PAR=M1PAR+VETT(1,COUNT) + M1PER=M1PER+VETT(2,COUNT) + CB=CB+1 + PRINT*,'A1 =',VETT(1,COUNT),'A2 =',VETT(2,COUNT) + END DO + M1PAR=M1PAR/CB + M1PER=M1PER/CB + + CB=0 + DEVPAR=0 + DO COUNT=1,CONFPRE + DEVPAR=DEVPAR+(VETT(1,COUNT)-M1PAR)**2 + END DO + DEVPAR=DEVPAR/CONFPRE + PRINT*,'Averaged new A1 =',M1PAR + PRINT*,'Averaged new A2 =',M1PER + PRINT*,'Standard deviation for new A1 =',SQRT(DEVPAR) + + CALL DECLAR('CHIAX','DP',' ',DUMMY2,M1PAR) + CALL DECLAR('CHIRH','DP',' ',DUMMY2,M1PER) + + END IF + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE FANTAERRORE (XPX,XPY,XPZ,XPMX,XPMY,XPMZ,NNAT,XOBS, + & XTOLPROT,METHOD,LOCFLG) +C ------------------------------------------------------------------ +C Calculates error on tensor parameters through Monte Carlo method +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + DOUBLE PRECISION APARA,APERP + COMMON /FRUN/ APERP,APARA + DOUBLE PRECISION ERRPARA(10000),ERRRH(10000) + DOUBLE PRECISION FEPERC + INTEGER FENUMC + COMMON /FERRORE/FEPERC,FENUMC + DOUBLE PRECISION FDEPERC + INTEGER FDENUMC + COMMON /FDERRORE/FDEPERC,FDENUMC + DOUBLE PRECISION XPX(1000),XPY(1000),XPZ(1000), + & XOBS(1000),XTOLPROT(1000) + CHARACTER NOMEFILE*132 + DOUBLE PRECISION XPMX(1000),XPMY(1000),XPMZ(1000),WPROT(1000) + INTEGER NNAT,METHOD,I,A,J,IA,AA,SEL + DOUBLE PRECISION RAND,MEDERRPARA,DEVERRPARA,DEVERRRH,MEDERRRH + LOGICAL LOCFLG + + IF (METHOD.EQ.0) THEN + FENUMC=FDENUMC + FEPERC=FDEPERC + END IF + NOMEFILE='MONTECARLO.OUT' + PRINT*,'FEPERC is:',FEPERC + + DO A=1,FENUMC + J=0 + DO I=1,NNAT + WPROT(I)=1 + END DO + AA=NINT(FEPERC*NNAT) + PRINT*,'AA is:',AA + IA=0 + DO WHILE (IA.LT.AA) + SEL=NINT(RAND()*DBLE(NNAT)) + PRINT*,'SEL IS:',SEL + IF (WPROT(SEL).EQ.1) THEN + WPROT(SEL)=0 + IA=IA+1 + END IF + END DO + CALL FANTALINEXE(XPX,XPY,XPZ,XPMX,XPMY,XPMZ,NNAT,XOBS,XTOLPROT, + & METHOD,NOMEFILE,WPROT,1,LOCFLG) + ERRPARA(A)=APARA + ERRRH(A)=APERP + PRINT*,'A1 =',APARA + PRINT*,'A2 =',APERP + END DO + + MEDERRPARA=0 + MEDERRRH=0 + DO A=1,FENUMC + MEDERRPARA=ERRPARA(A)+MEDERRPARA + MEDERRRH=ERRRH(A)+MEDERRRH + END DO + MEDERRPARA=MEDERRPARA/DBLE(FENUMC) + MEDERRRH=MEDERRRH/DBLE(FENUMC) + + DEVERRRH=0 + DEVERRPARA=0 + DO A=1,FENUMC + DEVERRPARA=DEVERRPARA+(ERRPARA(A)-MEDERRPARA)**2 + DEVERRRH=DEVERRRH+(ERRRH(A)-MEDERRRH)**2 + END DO + DEVERRPARA=DEVERRPARA/DBLE(FENUMC) + DEVERRRH=DEVERRRH/DBLE(FENUMC) + + PRINT*,'Standard deviation for A1 =',SQRT(DEVERRPARA) + PRINT*,'Standard deviation for A2 =',SQRT(DEVERRRH) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE EXPCS (EDI,WHICH) +C ------------------------------------------------------------------ +C Calls EXPCS2, which does the actual energy calculation +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "xdipo_pcs.fcm" + INCLUDE "heap.fcm" + + DOUBLE PRECISION EDI + CHARACTER*7 WHICH + + CALL EXPCS2 (EDI,HEAP(XPCSIPTR),HEAP(XPCSJPTR),HEAP(XPCSKPTR), + & HEAP(XPCSLPTR),HEAP(XPCSMPTR),HEAP(XPCSNPTR), + & HEAP(XPCSILST),HEAP(XPCSJLST),HEAP(XPCSKLST), + & HEAP(XPCSLLST),HEAP(XPCSMLST),HEAP(XPCSNLST), + & HEAP(XPCSOBSPTR),HEAP(XPCSERRPTR), + & HEAP(CALCXPCSPTR),WHICH) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE EXPCS2 (EDI,ATOMIPTR,ATOMJPTR,ATOMKPTR, + & ATOMLPTR,ATOMMPTR,ATOMNPTR, + & ATOMILST,ATOMJLST,ATOMKLST, + & ATOMLLST,ATOMMLST,ATOMNLST, + & XPCSOBS,XPCSERR,XPCSCALC,WHICH) +C ------------------------------------------------------------------ +C Calculates pseudocontact shift energies +C +C Energies are of the form: +C E = K*(DELTAPCS**2) +C where: +C K = FORCE CONSTANT +C DELTAPCS = CALCULATED PCS - OBSERVED PCS +C and pseudocontact shift function is defined as: +C PCS=A1*(3*COS(THETA)^2-1)+(3/2)*A2*(SIN(THETA)^2*COS(2*PHI))/R3 +C +C WHICH is a flag that switches between energy/force calculations +C and PCS calculations for violations +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "psf.fcm" + INCLUDE "coord.fcm" + INCLUDE "numbers.fcm" + INCLUDE "deriv.fcm" + INCLUDE "xdipo_pcs.fcm" + INCLUDE "consta.fcm" + + DOUBLE PRECISION XPCSE + COMMON /XPCSENERGY/XPCSE + INTEGER ATOMIPTR(*),ATOMJPTR(*),ATOMKPTR(*) + INTEGER ATOMLPTR(*),ATOMMPTR(*),ATOMNPTR(*) + INTEGER ATOMILST(*),ATOMJLST(*),ATOMKLST(*) + INTEGER ATOMLLST(*),ATOMMLST(*),ATOMNLST(*) + DOUBLE PRECISION XPCSOBS(*),XPCSERR(*) + DOUBLE PRECISION XPCSCALC(*),EDI + CHARACTER*7 WHICH +C ------------------------------------------------------------------ +C Local variables +C ------------------------------------------------------------------ + INTEGER COUNT,CLASS,MM,NN,II,I,J,K,L,NA,M,N + PARAMETER (NA=20) + INTEGER DMM,DNN,OUTFLAG(NA),MCOUNT,NCOUNT,PCOUNT + DOUBLE PRECISION XI,XJ,XK,XL,XM,XN,YI,YJ,YK,YL,YM,YN, + & ZI,ZJ,ZK,ZL,ZM,ZN, + & E,DF1(NA),DF2(NA),DF3(NA), + & PHI(NA),DPXI(NA),DPYI(NA),DPZI(NA), + & DPXJ(NA),DPYJ(NA),DPZJ(NA),DPXK(NA), + & DPYK(NA),DPZK(NA),DPXL(NA),DPYL(NA), + & DPZL(NA),DPXM(NA),DPYM(NA),DPZM(NA), + & DPXN(NA),DPYN(NA),DPZN(NA), + & THETA(NA),DTXI(NA),DTYI(NA),DTZI(NA), + & DTXJ(NA),DTYJ(NA),DTZJ(NA),DTXK(NA), + & DTYK(NA),DTZK(NA),DTXL(NA),DTYL(NA), + & DTZL(NA),DTXM(NA),DTYM(NA),DTZM(NA), + & DTXN(NA),DTYN(NA),DTZN(NA),DRXM(NA), + & DRYM(NA),DRZM(NA),DRXN(NA),DRYN(NA), + & DRZN(NA),DIST, + & O1,O2,O3,O4,O5,O6,O7,O8,O9,O10, + & O11,O12,O13,O14,O15,O16,O17,O18,O19, + & O20,O21,O22,O23,O24,O25,O26,O27,O28, + & O29,O30,O31,O32,O33,O34,O35,O36,O37, + & O38,O39,O40,O41,O42,O43,O44,O45,O46, + & O47,O48,O49,O50,O51,O52,O53,O54,O55, + & O56,O57,O58,O59,O60,O61,O62,O63,O64, + & O65,O66,O67,O68,O69,O70,O71,O72,O73, + & O74,O75,O76,O77,O78,O79,O80,O81,O82 + DOUBLE PRECISION OBSXPCS,ERRXPCS,K1, + & COEF1,COEF2,A,B,DELTA, + & DT,DP,DELTAXPCS, + & DD,DR,CALCXPCS +C ------------------------------------------------------------------ +C Zero out partial energy +C ------------------------------------------------------------------ + EDI=ZERO + + CLASS=1 + K1=XPCSFORCES(1) + COEF1=XPCSCOEF1(1) + COEF2=XPCSCOEF2(1) + + COUNT=0 + DO WHILE (COUNT.LT.PCSNUM) + COUNT=COUNT+1 +C ------------------------------------------------------------------ +C Reset individual E to zero +C ------------------------------------------------------------------ + E=0 + + DO WHILE ((XPCSASSNDX(CLASS).LT.COUNT).OR. + & (XPCSASSNDX(CLASS).EQ.0)) + CLASS=CLASS+1 + END DO + + IF (XPCSASSNDX(CLASS).EQ.XPCSASSNDX(CLASS-1)) THEN + COUNT=COUNT-1 + END IF + + K1=XPCSFORCES(CLASS) + COEF1=XPCSCOEF1(CLASS) + COEF2=XPCSCOEF2(CLASS) +C ------------------------------------------------------------------ +C Note there should only be one atom for I, J, K, and L +C ------------------------------------------------------------------ + I=ATOMIPTR(COUNT)+1 + J=ATOMJPTR(COUNT)+1 + K=ATOMKPTR(COUNT)+1 + L=ATOMLPTR(COUNT)+1 + + XI=X(ATOMILST(I)) + XJ=X(ATOMJLST(J)) + XK=X(ATOMKLST(K)) + XL=X(ATOMLLST(L)) + YI=Y(ATOMILST(I)) + YJ=Y(ATOMJLST(J)) + YK=Y(ATOMKLST(K)) + YL=Y(ATOMLLST(L)) + ZI=Z(ATOMILST(I)) + ZJ=Z(ATOMJLST(J)) + ZK=Z(ATOMKLST(K)) + ZL=Z(ATOMLLST(L)) + + OBSXPCS=XPCSOBS(COUNT) + ERRXPCS=XPCSERR(COUNT) +C ------------------------------------------------------------------ +C Initialize calculated pseudocontact shift and counter +C ------------------------------------------------------------------ + CALCXPCS=0 + II=0 + MCOUNT=0 + NCOUNT=0 +C ------------------------------------------------------------------ +C Check for correct permutations of paired atoms to get the metal- +C nucleus vector. This depends solely on the order of the assigned +C atoms. OUTFLAG=1 indicates that the permutation is not allowed. +C ------------------------------------------------------------------ + DMM=ATOMMPTR(COUNT+1)-ATOMMPTR(COUNT) + DNN=ATOMNPTR(COUNT+1)-ATOMNPTR(COUNT) + DO MM=ATOMMPTR(COUNT)+1,ATOMMPTR(COUNT+1) + MCOUNT=MCOUNT+1 + NCOUNT=0 + DO NN=ATOMNPTR(COUNT)+1,ATOMNPTR(COUNT+1) + NCOUNT=NCOUNT+1 + II=II+1 + IF ((DMM.GT.1).AND.(DNN.GT.1)) THEN + IF (MCOUNT.EQ.NCOUNT) THEN + OUTFLAG(II)=0 + ELSE + OUTFLAG(II)=1 + ENDIF + ELSE + OUTFLAG(II)=0 + END IF + END DO + END DO + + II=0 + PCOUNT=0 + DO MM=ATOMMPTR(COUNT)+1,ATOMMPTR(COUNT+1) + DO NN=ATOMNPTR(COUNT)+1,ATOMNPTR(COUNT+1) + II=II+1 + IF (OUTFLAG(II).NE.1) THEN + PCOUNT=PCOUNT+1 + XM=X(ATOMMLST(MM)) + XN=X(ATOMNLST(NN)) + YM=Y(ATOMMLST(MM)) + YN=Y(ATOMNLST(NN)) + ZM=Z(ATOMMLST(MM)) + ZN=Z(ATOMNLST(NN)) +C ------------------------------------------------------------------ +C Calculate THETA, PHI and derivatives with respect to X, Y, Z of +C atoms I, J, K, L, M, N. This was done with Mathematica by Nico! +C ------------------------------------------------------------------ + O1=-XI+XJ + O2=O1**2 + O3=-YI+YJ + O4=O3**2 + O5=-ZI+ZJ + O6=O5**2 + O7=O2+O4+O6 + O8=SQRT(O7) + O9=1/O8 + O10=-XM+XN + O11=O1*O10 + O12=-YM+YN + O13=O12*O3 + O14=-ZM+ZN + O15=O14*O5 + O16=O11+O13+O15 + O17=O10**2 + O18=O12**2 + O19=O14**2 + O20=O17+O18+O19 + O21=SQRT(O20) + O22=1/O21 + O23=1/O7 + O24=O16**2 + O25=1/O20 + O26=O23*O24*O25 + O27=SQRT(1.D0-O26) + O28=1/O27 + O29=XM-XN + O30=O7**(3.D0/2.D0) + O31=1/O30 + O32=O1*O16*O22*O31 + O33=YM-YN + O34=O16*O22*O3*O31 + O35=ZM-ZN + O36=O16*O22*O31*O5 + O37=O20**(3.D0/2.D0) + O38=1/O37 + O39=O10*O16*O38*O9 + O40=O12*O16*O38*O9 + O41=O14*O16*O38*O9 + O42=-XI+XK + O43=O42**2 + O44=-YI+YK + O45=O44**2 + O46=-ZI+ZK + O47=O46**2 + O48=O43+O45+O47 + O49=SQRT(O48) + O50=-XI+XL + O51=O50**2 + O52=-YI+YL + O53=O52**2 + O54=-ZI+ZL + O55=O54**2 + O56=O51+O53+O55 + O57=SQRT(O56) + O58=1/O57 + O59=O10*O42 + O60=O12*O44 + O61=O14*O46 + O62=O59+O60+O61 + O63=1/O62 + O64=O10*O50 + O65=O12*O52 + O66=O14*O54 + O67=O64+O65+O66 + O68=O62**2 + O69=1/O68 + O70=O56**(3.D0/2.D0) + O71=1/O70 + O72=O49*O50*O63*O67*O71 + O73=1/O49 + O74=O42*O58*O63*O67*O73 + O75=1/O56 + O76=O67**2 + O77=O48*O69*O75*O76 + O78=1/(1.D0+O77) + O79=O49*O52*O63*O67*O71 + O80=O44*O58*O63*O67*O73 + O81=O49*O54*O63*O67*O71 + O82=O46*O58*O63*O67*O73 + + THETA(II)=ACOS(O16*O22*O9) + + DTXI(II)=-(O28*(O32+O22*O29*O9)) + DTYI(II)=-(O28*(O34+O22*O33*O9)) + DTZI(II)=-(O28*(O36+O22*O35*O9)) + DTXJ(II)=-(O28*(-O32+O10*O22*O9)) + DTYJ(II)=-(O28*(-O34+O12*O22*O9)) + DTZJ(II)=-(O28*(-O36+O14*O22*O9)) + DTXK(II)=0 + DTYK(II)=0 + DTZK(II)=0 + DTXL(II)=0 + DTYL(II)=0 + DTZL(II)=0 + DTXM(II)=-(O28*(O39+O22*O9*(XI-XJ))) + DTYM(II)=-(O28*(O40+O22*O9*(YI-YJ))) + DTZM(II)=-(O28*(O41+O22*O9*(ZI-ZJ))) + DTXN(II)=-(O28*(-O39+O1*O22*O9)) + DTYN(II)=-(O28*(-O40+O22*O3*O9)) + DTZN(II)=-(O28*(-O41+O22*O5*O9)) + + PHI(II)=ATAN(O49*O58*O63*O67) + + DPXI(II)=O78*(O29*O49*O58*O63-O29*O49*O58*O67*O69+ + & O72-O74) + DPYI(II)=O78*(O33*O49*O58*O63-O33*O49*O58*O67*O69+ + & O79-O80) + DPZI(II)=O78*(O35*O49*O58*O63-O35*O49*O58*O67*O69+ + & O81-O82) + DPXJ(II)=0 + DPYJ(II)=0 + DPZJ(II)=0 + DPXK(II)=(-(O10*O49*O58*O67*O69)+O74)*O78 + DPYK(II)=O78*(-(O12*O49*O58*O67*O69)+O80) + DPZK(II)=O78*(-(O14*O49*O58*O67*O69)+O82) + DPXL(II)=(O10*O49*O58*O63-O72)*O78 + DPYL(II)=O78*(O12*O49*O58*O63-O79) + DPZL(II)=O78*(O14*O49*O58*O63-O81) + DPXM(II)=O78*(-(O49*O58*O67*O69*(XI-XK))+ + & O49*O58*O63*(XI-XL)) + DPYM(II)=O78*(-(O49*O58*O67*O69*(YI-YK))+ + & O49*O58*O63*(YI-YL)) + DPZM(II)=O78*(-(O49*O58*O67*O69*(ZI-ZK))+ + & O49*O58*O63*(ZI-ZL)) + DPXN(II)=(O49*O50*O58*O63-O42*O49*O58*O67*O69)*O78 + DPYN(II)=(O49*O52*O58*O63-O44*O49*O58*O67*O69)*O78 + DPZN(II)=(O49*O54*O58*O63-O46*O49*O58*O67*O69)*O78 +C ------------------------------------------------------------------ +C Calculate the distance between the two atoms M, N (metal-nucleus) +C and calculate the derivatives as well +C ------------------------------------------------------------------ + O1=-XM+XN + O2=O1**2 + O3=-YM+YN + O4=O3**2 + O5=-ZM+ZN + O6=O5**2 + O7=O2+O4+O6 + O8=SQRT(O7) + O9=1/O8 + O10=O1*O9 + O11=O3*O9 + O12=O5*O9 + DIST=O8 + DRXM(II)=-O10 + DRYM(II)=-O11 + DRZM(II)=-O12 + DRXN(II)=O10 + DRYN(II)=O11 + DRZN(II)=O12 +C ------------------------------------------------------------------ +C Calculate the pseudocontact shift +C ------------------------------------------------------------------ + A=COEF1 + B=COEF2 + CALCXPCS=(A*(3.*COS(THETA(II))*COS(THETA(II))-1.)+ + & B*(3./2.)*(SIN(THETA(II))*SIN(THETA(II))* + & COS(2.*PHI(II))))/(DIST**3.0) + END IF + END DO + END DO + + IF (WHICH.EQ.'ANALYZE') THEN + XPCSCALC(COUNT)=CALCXPCS + END IF + + DELTA=(CALCXPCS-OBSXPCS) +C ------------------------------------------------------------------ +C Adjust the deviation based on the error range +C ------------------------------------------------------------------ + IF ((DELTA.LT.0.000).AND.(ABS(DELTA).GT.ERRXPCS)) THEN + DELTAXPCS=DELTA+ERRXPCS + ELSE IF ((DELTA.GT.0.000).AND.(DELTA.GT.ERRXPCS)) THEN + DELTAXPCS=DELTA-ERRXPCS + ELSE + DELTAXPCS=0.0 + END IF + DD=DELTAXPCS + + II = 0 + PCOUNT = 0 + DO MM=ATOMMPTR(COUNT)+1,ATOMMPTR(COUNT+1) + DO NN=ATOMNPTR(COUNT)+1,ATOMNPTR(COUNT+1) + II = II+1 + IF (OUTFLAG(II).NE.1) THEN + PCOUNT=PCOUNT + 1 +C ------------------------------------------------------------------ +C Taking care of derivative +C ------------------------------------------------------------------ + DT=(1.0/(DIST**3))*(A*(-6.)*SIN(THETA(II))* + & COS(THETA(II))+B*3.*SIN(THETA(II))* + & COS(THETA(II))*COS(2.*PHI(II))) + DP=(1.0/DIST**3)*(B*(-3.)*SIN(THETA(II))* + & SIN(THETA(II))*SIN(2.*PHI(II))) + DR=(-3.0/(DIST**4))*(A*(3.*COS(THETA(II))* + & COS(THETA(II))-1.)+B*(3./2.)* + & (SIN(THETA(II))*SIN(THETA(II))* + & COS(2.*PHI(II)))) + IF (DD.EQ.0.0) THEN + E=E+0.0 + DF1(II)=0.0 + DF2(II)=0.0 + DF3(II)=0.0 + ELSE + E=E+K1*(DELTAXPCS**2) + DF1(II)=2*K1*DELTAXPCS*DT + DF2(II)=2*K1*DELTAXPCS*DP + DF3(II)=2*K1*DELTAXPCS*DR + END IF + END IF + END DO + END DO +C ------------------------------------------------------------------ +C Accumulate energy +C ------------------------------------------------------------------ + EDI=EDI+E +C ------------------------------------------------------------------ +C Now update forces if in energy/force mode +C ------------------------------------------------------------------ + IF (WHICH.NE.'ANALYZE') THEN + II = 0 + DO MM=ATOMMPTR(COUNT)+1,ATOMMPTR(COUNT+1) + DO NN=ATOMNPTR(COUNT)+1,ATOMNPTR(COUNT+1) + II=II+1 + IF (OUTFLAG(II).NE.1) THEN + DX(ATOMILST(I))=DX(ATOMILST(I))+ + & DF1(II)*DTXI(II)+DF2(II)*DPXI(II) + DY(ATOMILST(I))=DY(ATOMILST(I))+ + & DF1(II)*DTYI(II)+DF2(II)*DPYI(II) + DZ(ATOMILST(I))=DZ(ATOMILST(I))+ + & DF1(II)*DTZI(II)+DF2(II)*DPZI(II) + DX(ATOMJLST(J))=DX(ATOMJLST(J))+ + & DF1(II)*DTXJ(II)+DF2(II)*DPXJ(II) + DY(ATOMJLST(J))=DY(ATOMJLST(J))+ + & DF1(II)*DTYJ(II)+DF2(II)*DPYJ(II) + DZ(ATOMJLST(J))=DZ(ATOMJLST(J))+ + & DF1(II)*DTZJ(II)+DF2(II)*DPZJ(II) + DX(ATOMKLST(K))=DX(ATOMKLST(K))+ + & DF1(II)*DTXK(II)+DF2(II)*DPXK(II) + DY(ATOMKLST(K))=DY(ATOMKLST(K))+ + & DF1(II)*DTYK(II)+DF2(II)*DPYK(II) + DZ(ATOMKLST(K))=DZ(ATOMKLST(K))+ + & DF1(II)*DTZK(II)+DF2(II)*DPZK(II) + DX(ATOMLLST(L))=DX(ATOMLLST(L))+ + & DF1(II)*DTXL(II)+DF2(II)*DPXL(II) + DY(ATOMLLST(L))=DY(ATOMLLST(L))+ + & DF1(II)*DTYL(II)+DF2(II)*DPYL(II) + DZ(ATOMLLST(L))=DZ(ATOMLLST(L))+ + & DF1(II)*DTZL(II)+DF2(II)*DPZL(II) + DX(ATOMMLST(MM))=DX(ATOMMLST(MM))+ + & DF1(II)*DTXM(II)+DF2(II)* + & DPXM(II)+DF3(II)*DRXM(II) + DY(ATOMMLST(MM))=DY(ATOMMLST(MM))+ + & DF1(II)*DTYM(II)+DF2(II)* + & DPYM(II)+DF3(II)*DRYM(II) + DZ(ATOMMLST(MM))=DZ(ATOMMLST(MM))+ + & DF1(II)*DTZM(II)+DF2(II)* + & DPZM(II)+DF3(II)*DRZM(II) + DX(ATOMNLST(NN))=DX(ATOMNLST(NN))+ + & DF1(II)*DTXN(II)+DF2(II)* + & DPXN(II)+DF3(II)*DRXN(II) + DY(ATOMNLST(NN))=DY(ATOMNLST(NN))+ + & DF1(II)*DTYN(II)+DF2(II)* + & DPYN(II)+DF3(II)*DRYN(II) + DZ(ATOMNLST(NN))=DZ(ATOMNLST(NN))+ + & DF1(II)*DTZN(II)+DF2(II)* + & DPZN(II)+DF3(II)*DRZN(II) + + END IF + END DO + END DO + END IF + END DO + + XPCSE=EDI + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE READXPCS +C ------------------------------------------------------------------ +C Reads in pseudocontact shift information +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "comand.fcm" + INCLUDE "xdipo_pcs.fcm" + INCLUDE "funct.fcm" + INCLUDE "psf.fcm" + INCLUDE "heap.fcm" + INCLUDE "numbers.fcm" + + LOGICAL FSWITCH,FSAVE,FERR + COMMON /FPARAM/FSWITCH,FSAVE,FERR + INTEGER FVIOLS + COMMON /VIOL/FVIOLS + DOUBLE PRECISION FEPERC + INTEGER FENUMC + COMMON /FERRORE/FEPERC,FENUMC + + INTEGER COUNT,SPTR,OLDCLASS,OLDMAXXPCS,FNCLASS + DOUBLE PRECISION K1,CUTOFF,COEF1,COEF2 + DOUBLE PRECISION FMEDPERC + CHARACTER*4 THENAME + CHARACTER*132 NOMEFILE + INTEGER TOLSWI + + NOMEFILE='XDIPO_PCS.OUT' + + SPTR=ALLHP(INTEG4(NATOM)) + CALL PUSEND('XDIPO_PCS>') +862 CONTINUE + CALL NEXTWD('XDIPO_PCS>') + CALL MISCOM('XDIPO_PCS>',USED) + IF (.NOT.USED) THEN +C ------------------------------------------------------------------ +C Documentation +C ------------------------------------------------------------------ + IF (WD(1:4).EQ.'HELP') THEN + WRITE(DUNIT,'(10X,A)') + & ' XPCS {} END ', + & ' :== ', + & ' ASSIgn ', + & ' ', + & ' * Restraint: Metal Z X Y Nucleus PCS Err *', + & ' CLASsification ', + & ' * Starts a new class *', + & ' TOLL <1|0>', + & ' * Switch tolerance in FANTALIN (1=Yes, 0=No)', + & ' COEFficient ', + & ' * Tensor parameters: A1 A2 *', + & ' FORCeconstant ', + & ' * Force constant for the current class *', + & ' NREStraints ', + & ' * Number of slots for restraints to allocate *', + & ' PRINt THREshold ', + & ' * Prints violations larger than the threshold *', + & ' RESEt', + & ' * Erases the restraint table, but keeps NRES *', + & ' SAVE', + & ' * Filename to save PCS values *', + & ' FMED', + & ' * Calculates average tensor parameters *', + & ' ERRON', + & ' * Switches Monte Carlo error evaluation ON *', + & ' ERROFF', + & ' * Switches Monte Carlo error evaluation OFF *', + & ' FON', + & ' * Switches tensor auto-update mode ON *', + & ' FOFF', + & ' * Switches tensor auto-update mode OFF *', + & ' SON', + & ' * Switches saving mode ON *', + & ' SOFF', + & ' * Switches saving mode OFF *', + & ' FRUN', + & ' * Runs FANTALIN *' +C ------------------------------------------------------------------ +C About FANTALIN +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'TOLL') THEN + CALL NEXTI('Tolerance setting in FANTALIN (1/0) =',TOLSWI) + IF (TOLSWI.EQ.1) THEN + PCSFLG=.TRUE. + WRITE (DUNIT,'(A)') 'Tolerance ON in FANTALIN.' + ELSE IF (TOLSWI.EQ.0) THEN + PCSFLG=.FALSE. + WRITE (DUNIT,'(A)') 'Tolerance OFF in FANTALIN.' + ELSE + WRITE (DUNIT,'(A)') 'Unknown switch. Using default.' + END IF + ELSE IF (WD(1:4).EQ.'SAVE') THEN + CALL NEXTFI('FILENAME =',NOMEFILE) + PRINT*,'Name of the file to save PCS values :',NOMEFILE + ELSE IF (WD(1:3).EQ.'FME') THEN + CALL NEXTF('Percentage for average (0/100) =',FMEDPERC) + CALL NEXTI('Average on class number =',FNCLASS) + IF (FNCLASS.GT.NCLASSES) THEN + PRINT*,'%FRUN-ERR: This class does not exist...' + ELSE + CALL FANTAMEDIA(FNCLASS,FMEDPERC) + END IF + ELSE IF (WD(1:5).EQ.'ERRON') THEN + CALL NEXTI('Number of cycles =',FENUMC) + CALL NEXTF('Percentage to discard (0/100) =',FEPERC) + FEPERC=FEPERC/100.0 + FERR=1 + PRINT*,'Monte Carlo error evaluation ON' + ELSE IF (WD(1:6).EQ.'ERROFF') THEN + FERR=0 + PRINT*,'Monte Carlo error evaluation OFF' + ELSE IF (WD(1:3).EQ.'FON') THEN + FSWITCH=1 + PRINT*,'Tensor auto-update mode ON' + ELSE IF (WD(1:4).EQ.'FOFF') THEN + FSWITCH=0 + PRINT*,'Tensor auto-update mode OFF' + ELSE IF (WD(1:3).EQ.'SON') THEN + FSAVE=1 + PRINT*,'Saving mode ON' + ELSE IF (WD(1:4).EQ.'SOFF') THEN + FSAVE=0 + PRINT*,'Saving mode OFF' + ELSE IF (WD(1:4).EQ.'FRUN') THEN + CALL NEXTI('FANTALIN on class number =',FNCLASS) + IF (FNCLASS.GT.NCLASSES.OR.FNCLASS.EQ.0) THEN + PRINT*,'%FRUN-ERR: This class does not exist...' + ELSE + CALL FANTALIN(HEAP(XPCSIPTR),HEAP(XPCSJPTR), + & HEAP(XPCSKPTR),HEAP(XPCSLPTR), + & HEAP(XPCSMPTR),HEAP(XPCSNPTR), + & HEAP(XPCSILST),HEAP(XPCSJLST), + & HEAP(XPCSKLST),HEAP(XPCSLLST), + & HEAP(XPCSMLST),HEAP(XPCSNLST), + & HEAP(XPCSOBSPTR),HEAP(XPCSERRPTR), + & FNCLASS,NOMEFILE ) + END IF +C ------------------------------------------------------------------ +C Get class name. Determine if it's an already-defined class. +C Insert a new class if it's not. +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'CLAS') THEN + OLDCLASS=CURCLASS + CALL NEXTA4('Class name =',THENAME) + MODE=NEW + DO COUNT=1,NCLASSES + IF (XPCSCLASSNAMES(COUNT).EQ.THENAME) THEN + MODE=UPDATE + CURCLASS=COUNT + END IF + END DO + IF (MODE.EQ.NEW) THEN +C ------------------------------------------------------------------ +C Make sure you can't add more than the maximum number of classes +C ------------------------------------------------------------------ + IF (OLDCLASS.EQ.MAXXPCSCLASSES) THEN + CALL DSPERR('XDIPO_PCS','Too many classes.') + CALL DSPERR('XDIPO_PCS', + & 'Increase MAXXPCSCLASSES and recompile.') + CALL WRNDIE(-5,'READXPCS', + & 'Too many anisotropy classes.') + END IF + NCLASSES=NCLASSES+1 + CURCLASS=NCLASSES + XPCSCLASSNAMES(CURCLASS)=THENAME + +C ------------------------------------------------------------------ +C If this isn't the first class, close off the old class +C ------------------------------------------------------------------ + IF (NCLASSES.GT.1) THEN + XPCSASSNDX(OLDCLASS)=PCSNUM + END IF + END IF +C ------------------------------------------------------------------ +C Set force constant for current class +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'FORC') THEN +C ------------------------------------------------------------------ +C Start a default class if there isn't one defined +C ------------------------------------------------------------------ + IF (CURCLASS.EQ.0) THEN + NCLASSES=1 + CURCLASS=1 + END IF + CALL NEXTF('Force constant =',K1) + WRITE(DUNIT,'(A,A,A,F8.3)') + & 'Setting force constant for class ', + & XPCSCLASSNAMES(CURCLASS),' to ',K1 + XPCSFORCES(CURCLASS)=K1 +C ------------------------------------------------------------------ +C Set coefficient constant for current class +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'COEF') THEN + CALL NEXTF('Coefficient A1 =',COEF1) + CALL NEXTF('Coefficient A2 =',COEF2) +C ------------------------------------------------------------------ +C Start a default class if there isn't one defined +C ------------------------------------------------------------------ + IF (CURCLASS.EQ.0) THEN + NCLASSES=1 + CURCLASS=1 + END IF + WRITE(DUNIT,'(A,A,A,F8.3,F8.3)') + & 'Setting coefficients for class ', + & XPCSCLASSNAMES(CURCLASS),' to ',COEF1,COEF2 + XPCSCOEF1(CURCLASS)=COEF1 + XPCSCOEF2(CURCLASS)=COEF2 +C ------------------------------------------------------------------ +C Reset pseudocontact shifts database +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'RESE') THEN + CALL XPCSDEFAULTS + CALL ALLOCXPCS(0,MAXXPCS) +C ------------------------------------------------------------------ +C Change number of assignment slots +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'NRES') THEN + OLDMAXXPCS=MAXXPCS + CALL NEXTI('Number of slots =',MAXXPCS) + CALL ALLOCXPCS(OLDMAXXPCS,MAXXPCS) + WRITE(DUNIT,'(A,I8,A)') + & 'XDIPO_PCS: Allocating space for',MAXXPCS, + & 'number of restraints.' +C ------------------------------------------------------------------ +C Read in an assignment +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'ASSI') THEN +C ------------------------------------------------------------------ +C Make sure you can't add more assignments than you have slots for +C ------------------------------------------------------------------ + IF (XPCSMX.EQ.MAXXPCS) THEN + CALL DSPERR('XDIPO_PCS','Too many assignments.') + CALL DSPERR('XDIPO_PCS', + & 'Increasing NREStraints by 100.') + OLDMAXXPCS=MAXXPCS + MAXXPCS=MAXXPCS+100 + CALL ALLOCXPCS(OLDMAXXPCS,MAXXPCS) + END IF +C ------------------------------------------------------------------ +C If there isn't a class specified, start a default class +C ------------------------------------------------------------------ + IF (CURCLASS.EQ.0) THEN + NCLASSES=1 + CURCLASS=1 + END IF + CALL READXPCS2(HEAP(XPCSIPTR),HEAP(XPCSJPTR), + & HEAP(XPCSKPTR),HEAP(XPCSLPTR), + & HEAP(XPCSMPTR),HEAP(XPCSNPTR), + & HEAP(XPCSILST),HEAP(XPCSJLST), + & HEAP(XPCSKLST),HEAP(XPCSLLST), + & HEAP(XPCSMLST),HEAP(XPCSNLST), + & HEAP(XPCSOBSPTR),HEAP(XPCSERRPTR), + & HEAP(SPTR)) +C ------------------------------------------------------------------ +C Print violations +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'PRIN') THEN + CALL NEXTWD('PRINT>') + IF (WD(1:4).NE.'THRE') THEN + CALL DSPERR('XDIPO_PCS', + & 'Print expects THREshold parameter.') + ELSE + CALL NEXTF('THREshold =',CUTOFF) + IF (CUTOFF.LT.ZERO) THEN + CALL DSPERR('XDIPO_PCS', + & 'Cutoff must be positive.') + CUTOFF=ABS(CUTOFF) + END IF + CALL NEXTA4('ALL OR CLASs>',THENAME) + IF (THENAME(1:3).EQ.'ALL') THEN + DO COUNT=1,NCLASSES + PRINTCLASS(COUNT)=.TRUE. + END DO + ELSE IF (THENAME(1:4).EQ.'CLAS') THEN + CALL NEXTA4('CLASS NAME =',THENAME) + DO COUNT=1,NCLASSES + IF (XPCSCLASSNAMES(COUNT).EQ.THENAME) THEN + PRINTCLASS(COUNT)=.TRUE. + ELSE + PRINTCLASS(COUNT)=.FALSE. + END IF + END DO + ELSE + DO COUNT=1,NCLASSES + PRINTCLASS(COUNT)=.TRUE. + END DO + END IF + CALL PRINTXPCS(CUTOFF,HEAP(CALCXPCSPTR), + & HEAP(XPCSOBSPTR),HEAP(XPCSERRPTR), + & HEAP(XPCSIPTR),HEAP(XPCSJPTR), + & HEAP(XPCSKPTR),HEAP(XPCSLPTR), + & HEAP(XPCSMPTR),HEAP(XPCSNPTR), + & HEAP(XPCSILST),HEAP(XPCSJLST), + & HEAP(XPCSKLST),HEAP(XPCSLLST), + & HEAP(XPCSMLST),HEAP(XPCSNLST)) + END IF +C ------------------------------------------------------------------ +C Check for END statement +C ------------------------------------------------------------------ + ELSE + CALL CHKEND('XDIPO_PCS>',DONE) + END IF + END IF + IF (.NOT.DONE) GOTO 862 + DONE=.FALSE. + CALL FREHP(SPTR,INTEG4(NATOM)) + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE ALLOCXPCS (OLDSIZE,NEWSIZE) +C ------------------------------------------------------------------ +C Reset pseudocontact shift arrays to hold size entries +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "funct.fcm" + INCLUDE "xdipo_pcs.fcm" + + INTEGER OLDSIZE,NEWSIZE + + IF (OLDSIZE.NE.0) THEN + CALL FREHP(XPCSIPTR,INTEG4(OLDSIZE)) + CALL FREHP(XPCSJPTR,INTEG4(OLDSIZE)) + CALL FREHP(XPCSKPTR,INTEG4(OLDSIZE)) + CALL FREHP(XPCSLPTR,INTEG4(OLDSIZE)) + CALL FREHP(XPCSMPTR,INTEG4(OLDSIZE)) + CALL FREHP(XPCSNPTR,INTEG4(OLDSIZE)) + + CALL FREHP(XPCSILST,INTEG4(OLDSIZE)) + CALL FREHP(XPCSJLST,INTEG4(OLDSIZE)) + CALL FREHP(XPCSKLST,INTEG4(OLDSIZE)) + CALL FREHP(XPCSLLST,INTEG4(OLDSIZE)) + CALL FREHP(XPCSMLST,INTEG4(OLDSIZE)) + CALL FREHP(XPCSNLST,INTEG4(OLDSIZE)) + + CALL FREHP(XPCSOBSPTR,IREAL8(OLDSIZE)) + CALL FREHP(XPCSERRPTR,IREAL8(OLDSIZE)) + CALL FREHP(CALCXPCSPTR,IREAL8(OLDSIZE)) + END IF + XPCSIPTR=ALLHP(INTEG4(NEWSIZE)) + XPCSJPTR=ALLHP(INTEG4(NEWSIZE)) + XPCSKPTR=ALLHP(INTEG4(NEWSIZE)) + XPCSLPTR=ALLHP(INTEG4(NEWSIZE)) + XPCSMPTR=ALLHP(INTEG4(NEWSIZE)) + XPCSNPTR=ALLHP(INTEG4(NEWSIZE)) + + XPCSILST=ALLHP(INTEG4(NEWSIZE)) + XPCSJLST=ALLHP(INTEG4(NEWSIZE)) + XPCSKLST=ALLHP(INTEG4(NEWSIZE)) + XPCSLLST=ALLHP(INTEG4(NEWSIZE)) + XPCSMLST=ALLHP(INTEG4(NEWSIZE)) + XPCSNLST=ALLHP(INTEG4(NEWSIZE)) + + XPCSOBSPTR=ALLHP(IREAL8(NEWSIZE)) + XPCSERRPTR=ALLHP(IREAL8(NEWSIZE)) + CALCXPCSPTR=ALLHP(IREAL8(NEWSIZE)) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE XPCSDEFAULTS +C ------------------------------------------------------------------ +C Sets up defaults +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "xdipo_pcs.fcm" + + INTEGER COUNT + + MODE=NEW + MAXXPCS=200 + XPCSMX=200 + PCSNUM=0 + NCLASSES=0 + CURCLASS=0 + PCSFLG=.FALSE. + DO COUNT=1,MAXXPCSCLASSES + XPCSCLASSNAMES(COUNT)='DEFAULT' + XPCSASSNDX(COUNT)=0 + XPCSFORCES(COUNT)=5.0 + XPCSCOEF1(COUNT)=1.0 + XPCSCOEF2(COUNT)=1.0 + END DO + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE READXPCS2(ATOMIPTR,ATOMJPTR,ATOMKPTR,ATOMLPTR,ATOMMPTR, + & ATOMNPTR,ATOMILST,ATOMJLST,ATOMKLST,ATOMLLST, + & ATOMMLST,ATOMNLST,XPCSOBS,XPCSERR,SEL) +C ------------------------------------------------------------------ +C Reads actual pseudocontact shift assignments into arrays +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "psf.fcm" + INCLUDE "coord.fcm" + INCLUDE "xdipo_pcs.fcm" + + INTEGER ATOMIPTR(*),ATOMJPTR(*),ATOMKPTR(*),ATOMLPTR(*) + INTEGER ATOMMPTR(*),ATOMNPTR(*),ATOMILST(*) + INTEGER ATOMJLST(*),ATOMKLST(*),ATOMLLST(*),ATOMMLST(*) + INTEGER ATOMNLST(*),SEL(*) + INTEGER ITMP,JTMP,KTMP,LTMP,II + DOUBLE PRECISION XPCSOBS(*),XPCSERR(*) +C ------------------------------------------------------------------ +C Local variables +C ------------------------------------------------------------------ + INTEGER NSEL,INSERTPOS,COUNT,CURSTOP,OTHERSTOP,NFLAG + INTEGER I + DOUBLE PRECISION XPCSO,XPCSE +C ------------------------------------------------------------------ +C If we're in update mode, make a space for the new line +C ------------------------------------------------------------------ + NFLAG=0 + IF (MODE.EQ.UPDATE) THEN + DO COUNT=PCSNUM+1,XPCSASSNDX(CURCLASS)+1,-1 + ATOMIPTR(COUNT)=ATOMIPTR(COUNT-1) + ATOMJPTR(COUNT)=ATOMJPTR(COUNT-1) + ATOMKPTR(COUNT)=ATOMKPTR(COUNT-1) + ATOMLPTR(COUNT)=ATOMLPTR(COUNT-1) + ATOMMPTR(COUNT)=ATOMMPTR(COUNT-1) + ATOMNPTR(COUNT)=ATOMNPTR(COUNT-1) + XPCSOBS(COUNT)=XPCSOBS(COUNT-1) + XPCSERR(COUNT)=XPCSERR(COUNT-1) + END DO + CURSTOP=XPCSASSNDX(CURCLASS) + DO COUNT=1,NCLASSES + OTHERSTOP=XPCSASSNDX(COUNT) + IF (OTHERSTOP.GT.CURSTOP) THEN + XPCSASSNDX(COUNT)=OTHERSTOP+1 + END IF + END DO + XPCSASSNDX(CURCLASS)=CURSTOP+1 + INSERTPOS=CURSTOP + PCSNUM=PCSNUM+1 + ELSE + PCSNUM=PCSNUM+1 + INSERTPOS=PCSNUM + XPCSASSNDX(CURCLASS)=INSERTPOS + END IF +C ------------------------------------------------------------------ +C Reading in the atom selection in the restraint table +C ------------------------------------------------------------------ + CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.) + IF (NSEL.GT.1) THEN + CALL DSPERR('XDIPO_PCS', + & 'More than 1 atom in SEL for atom I. Using first') + END IF + IF (NSEL.EQ.0) THEN + CALL DSPERR('XDIPO_PCS','Error with atom selection') + NFLAG=1 + END IF + CALL MAKIND(SEL,NATOM,NSEL) + IF (INSERTPOS.EQ.1) ATOMIPTR(INSERTPOS)=0 + II=ATOMIPTR(INSERTPOS) + II=II+1 + IF (II.GT.XPCSMX) XPCSMX=II + ATOMILST(II)=SEL(1) + ATOMIPTR(INSERTPOS+1)=II +C ------------------------------------------------------------------ +C Both the atoms I and M are the metal center +C ------------------------------------------------------------------ + IF (INSERTPOS.EQ.1) ATOMMPTR(INSERTPOS)=0 + II=ATOMMPTR(INSERTPOS) + II=II+1 + IF (II.GT.XPCSMX) XPCSMX=II + ATOMMLST(II)=SEL(1) + ATOMMPTR(INSERTPOS+1)=II + + CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.) + IF (NSEL.GT.1) THEN + CALL DSPERR('XDIPO_PCS', + & 'More than 1 atom in SEL for atom J. Using first') + END IF + IF (NSEL.EQ.0) THEN + CALL DSPERR('XDIPO_PCS','Error with atom selection') + NFLAG=1 + END IF + CALL MAKIND(SEL,NATOM,NSEL) + IF (INSERTPOS.EQ.1) ATOMJPTR(INSERTPOS)=0 + II=ATOMJPTR(INSERTPOS) + II=II+1 + IF (II.GT.XPCSMX) XPCSMX=II + ATOMJLST(II)=SEL(1) + ATOMJPTR(INSERTPOS+1)=II + + CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.) + IF (NSEL.GT.1) THEN + CALL DSPERR('XDIPO_PCS', + & 'More than 1 atom in SEL for atom K. Using first') + END IF + IF (NSEL.EQ.0) THEN + CALL DSPERR('XDIPO_PCS','Error with atom selection') + NFLAG=1 + END IF + CALL MAKIND(SEL, NATOM, NSEL) + IF (INSERTPOS.EQ.1) ATOMKPTR(INSERTPOS)=0 + II=ATOMKPTR(INSERTPOS) + II=II+1 + IF (II.GT.XPCSMX) XPCSMX=II + ATOMKLST(II)=SEL(1) + ATOMKPTR(INSERTPOS+1)=II + + CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.) + IF (NSEL.GT.1) THEN + CALL DSPERR('XDIPO_PCS', + & 'More than 1 atom in SEL for atom L. Using first') + END IF + IF (NSEL.EQ.0) THEN + CALL DSPERR('XDIPO_PCS','Error with atom selection') + NFLAG=1 + END IF + CALL MAKIND(SEL,NATOM,NSEL) + IF (INSERTPOS.EQ.1) ATOMLPTR(INSERTPOS)=0 + II=ATOMLPTR(INSERTPOS) + II=II+1 + IF (II.GT.XPCSMX) XPCSMX=II + ATOMLLST(II)=SEL(1) + ATOMLPTR(INSERTPOS+1)=II +C ------------------------------------------------------------------ +C Next selection defines the nucleus whose PCS was measured +C ------------------------------------------------------------------ + CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.) + IF (NSEL.GT.1) THEN + CALL DSPERR('XDIPO_PCS', + & 'More than 1 atom in SEL for atom N. Using first') + END IF + IF (NSEL.EQ.0) THEN + CALL DSPERR('XDIPO_PCS','Error with atom selection') + NFLAG=1 + END IF + CALL MAKIND(SEL,NATOM,NSEL) + IF (INSERTPOS.EQ.1) ATOMNPTR(INSERTPOS)=0 + II=ATOMNPTR(INSERTPOS) + II=II+1 + IF (II.GT.XPCSMX) XPCSMX=II + ATOMNLST(II)=SEL(1) + ATOMNPTR(INSERTPOS+1)=II +C ------------------------------------------------------------------ +C Reading in the observed pseudocontact shift +C ------------------------------------------------------------------ + CALL NEXTF('Observed PCS =', XPCSO) + CALL NEXTF('Error in PCS =', XPCSE) + + XPCSOBS(INSERTPOS)=XPCSO + XPCSERR(INSERTPOS)=XPCSE +C ------------------------------------------------------------------ +C Check for error atom selection. If there is one then reset the +C counter for restraint +C ------------------------------------------------------------------ + IF (NFLAG.EQ.1) THEN + PCSNUM=PCSNUM-1 + END IF + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE XPCSINIT +C ------------------------------------------------------------------ +C Initializes pseudocontact shift stuff +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "xdipo_pcs.fcm" + + CALL XPCSDEFAULTS + CALL ALLOCXPCS(0,MAXXPCS) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE PRINTXPCS (CUTOFF,XPCSCALC,XPCSOBS,XPCSERR, + & ATOMIPTR,ATOMJPTR,ATOMKPTR, + & ATOMLPTR,ATOMMPTR,ATOMNPTR, + & ATOMILST,ATOMJLST,ATOMKLST, + & ATOMLLST,ATOMMLST,ATOMNLST) +C ------------------------------------------------------------------ +C Prints pseudocontact shifts with deviation larger than cutoff, +C calculates RMSD and puts it into $RESULT +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "xdipo_pcs.fcm" + INCLUDE "comand.fcm" + INCLUDE "psf.fcm" + INCLUDE "numbers.fcm" + + INTEGER FVIOLS + COMMON /VIOL/FVIOLS + DOUBLE PRECISION CUTOFF,XPCSCALC(*),XPCSOBS(*),XPCSERR(*) + INTEGER ATOMILST(*),ATOMJLST(*),ATOMKLST(*) + INTEGER ATOMLLST(*),ATOMMLST(*),ATOMNLST(*) + INTEGER ATOMIPTR(*),ATOMJPTR(*),ATOMKPTR(*) + INTEGER ATOMLPTR(*),ATOMMPTR(*),ATOMNPTR(*) +C ------------------------------------------------------------------ +C Local variables +C ------------------------------------------------------------------ + DOUBLE PRECISION CALCXPCS,OBSXPCS,DELTAXPCS,DELTA,DP + INTEGER COUNT,CLASS,I,J,K,L,M,N,NUM,II + DOUBLE PRECISION RMS,VIOLS,ERRXPCS,XPCSENERGY + DOUBLE COMPLEX DUMMY2 + LOGICAL PRINTTHISCLASS + + RMS=ZERO + VIOLS=ZERO + NUM=0 +C ------------------------------------------------------------------ +C Make sure that the array of calculated PCS is up to date +C ------------------------------------------------------------------ + CALL EXPCS(XPCSENERGY,'ANALYZE') + WRITE (PUNIT,'(A)') 'The following pseudocontact shifts have' + WRITE (PUNIT,'(A)') 'deviation larger than or' + WRITE (PUNIT,'(A)') 'equal to the cutoff:' +C ------------------------------------------------------------------ +C Write out first class heading +C ------------------------------------------------------------------ + CLASS=1 + PRINTTHISCLASS=PRINTCLASS(CLASS) + IF (PRINTTHISCLASS) THEN + WRITE (PUNIT,'(A,A)') 'Class ',XPCSCLASSNAMES(1) + END IF +C ------------------------------------------------------------------ +C For every pseudocontact shift entry... +C ------------------------------------------------------------------ + COUNT = 0 + DO WHILE (COUNT.LT.PCSNUM) + COUNT=COUNT+1 +C ------------------------------------------------------------------ +C Is this the start of a new class? +C ------------------------------------------------------------------ + IF (XPCSASSNDX(CLASS).LT.COUNT) THEN + CLASS=CLASS+1 + PRINTTHISCLASS=PRINTCLASS(CLASS) + IF (XPCSASSNDX(CLASS).EQ.XPCSASSNDX(CLASS-1)) THEN + COUNT=COUNT-1 + END IF + IF (PRINTTHISCLASS) THEN + WRITE (PUNIT,'(A,A)') 'Class ',XPCSCLASSNAMES(CLASS) + END IF + END IF +C ------------------------------------------------------------------ +C If this assignment is in a class to be printed +C and make sure there is an entry for that class +C ------------------------------------------------------------------ + IF ((PRINTTHISCLASS).AND. + & (XPCSASSNDX(CLASS).NE.XPCSASSNDX(CLASS-1))) THEN +C ------------------------------------------------------------------ +C Always update RMSD +C ------------------------------------------------------------------ + CALCXPCS=XPCSCALC(COUNT) + OBSXPCS=XPCSOBS(COUNT) + ERRXPCS=XPCSERR(COUNT) + DP=(CALCXPCS-OBSXPCS) + + IF ((DP.LT.0.000).AND.(ABS(DP).GT.ERRXPCS)) THEN + DELTAXPCS=DP+ERRXPCS + ELSE IF ((DP.GT.0.000).AND.(DP.GT.ERRXPCS)) THEN + DELTAXPCS=DP-ERRXPCS + ELSE + DELTAXPCS=0.0 + END IF + + RMS=RMS+DELTAXPCS**2 + NUM=NUM+1 +C ------------------------------------------------------------------ +C Print out deviations larger than cutoff +C and update number of violations +C ------------------------------------------------------------------ + IF (ABS(DELTAXPCS).GE.CUTOFF) THEN + I=ATOMILST(ATOMIPTR(COUNT)+1) + J=ATOMJLST(ATOMJPTR(COUNT)+1) + K=ATOMKLST(ATOMKPTR(COUNT)+1) + L=ATOMLLST(ATOMLPTR(COUNT)+1) + WRITE(PUNIT,'(A,A)') '===============================', + & '===============================' + WRITE(PUNIT,'(A)') ' Set-M-atoms' + DO II=ATOMMPTR(COUNT)+1,ATOMMPTR(COUNT+1) + M=ATOMMLST(II) + WRITE(PUNIT,'(9X,4(1X,A))') SEGID(M),RESID(M),RES(M), + & TYPE(M) + END DO + WRITE(PUNIT,'(A)') ' Set-N-atoms' + DO II=ATOMNPTR(COUNT)+1,ATOMNPTR(COUNT+1) + N=ATOMNLST(II) + WRITE(PUNIT,'(9X,4(1X,A))') SEGID(N),RESID(N),RES(N), + & TYPE(N) + END DO + + WRITE(PUNIT,'(2(2X,A,1X,F8.3))') + & 'Calc: ',CALCXPCS,'Obs: ',OBSXPCS + WRITE(PUNIT,'(2X,A,1X,F8.3,2X,A,1X,F8.3)') + & 'Error: ',ERRXPCS,'Delta: ',DELTAXPCS + VIOLS = VIOLS + ONE + END IF + END IF + END DO + + IF (NUM.GT.0) THEN + RMS=SQRT(RMS/NUM) + ELSE + RMS=0.0 + END IF + + FVIOLS=VIOLS + PRINT *,'FVIOLS =',FVIOLS + + CALL DECLAR('RESULT','DP',' ',DUMMY2,RMS) + CALL DECLAR('VIOLATIONS','DP',' ',DUMMY2,VIOLS) + + RETURN + END diff -u -N -r xplor-nih-2.9.2/source/xdipo_pcs.fcm PARAxplor-nih-2.9.2/source/xdipo_pcs.fcm --- xplor-nih-2.9.2/source/xdipo_pcs.fcm 1970-01-01 01:00:00.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/xdipo_pcs.fcm 2003-12-05 14:34:36.000000000 +0100 @@ -0,0 +1,61 @@ +C ------------------------------------------------------------------ +C Pseudocontact shift stuff +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + INTEGER MAXXPCSCLASSES + PARAMETER(MAXXPCSCLASSES=40) +C ------------------------------------------------------------------ +C Arrays that hold pseudocontact shift info +C XPCSASSNDX tells ending index of the PCS arrays for each class +C XPCSFORCES holds K1 for each class +C ------------------------------------------------------------------ + INTEGER XPCSASSNDX(MAXXPCSCLASSES) + REAL XPCSFORCES(MAXXPCSCLASSES) + REAL XPCSCOEF1(MAXXPCSCLASSES),XPCSCOEF2(MAXXPCSCLASSES) + CHARACTER*8 XPCSCLASSNAMES(MAXXPCSCLASSES) + LOGICAL PRINTCLASS(MAXXPCSCLASSES),PCSFLG +C ------------------------------------------------------------------ +C MAXXPCS is the number of slots set aside for PCS assignments +C PCSNUM is the total number of PCS entered +C ------------------------------------------------------------------ + INTEGER MAXXPCS,PCSNUM,NCLASSES,CURCLASS,XPCSMX +C ------------------------------------------------------------------ +C Pointers to arrays to hold atom numbers, observed PCS and errors +C ------------------------------------------------------------------ + INTEGER XPCSIPTR,XPCSJPTR,XPCSKPTR,XPCSLPTR,XPCSMPTR,XPCSNPTR, + & XPCSILST,XPCSJLST,XPCSKLST,XPCSLLST,XPCSMLST,XPCSNLST, + & XPCSOBSPTR,XPCSERRPTR,CALCXPCSPTR +C ------------------------------------------------------------------ +C Input modes +C ------------------------------------------------------------------ + INTEGER MODE,NEW,UPDATE + PARAMETER(NEW = 1) + PARAMETER(UPDATE = 2) +C ------------------------------------------------------------------ +C Parameters as set up in ETOR - Not used indeed +C ------------------------------------------------------------------ + DOUBLE PRECISION MCONST + PARAMETER(MCONST=0.0001D0) + DOUBLE PRECISION EPS + PARAMETER(EPS=0.1D0) +C ------------------------------------------------------------------ +C Common blocks +C ------------------------------------------------------------------ + COMMON /CXPCS/XPCSCLASSNAMES + COMMON /IXPCS1/XPCSASSNDX,MAXXPCS,PCSNUM,CURCLASS,NCLASSES,XPCSMX + COMMON /IXPCS2/XPCSIPTR,XPCSJPTR,XPCSKPTR, + & XPCSLPTR,XPCSMPTR,XPCSNPTR, + & XPCSILST,XPCSJLST,XPCSKLST, + & XPCSLLST,XPCSMLST,XPCSNLST, + & XPCSOBSPTR,XPCSERRPTR,CALCXPCSPTR,MODE + COMMON /RXPCS/XPCSFORCES,XPCSCOEF1,XPCSCOEF2 + COMMON /LXPCS/PRINTCLASS + COMMON /WXPCS/PCSFLG + + SAVE /CXPCS/ + SAVE /IXPCS1/ + SAVE /IXPCS2/ + SAVE /RXPCS/ + SAVE /LXPCS/ + SAVE /WXPCS/ diff -u -N -r xplor-nih-2.9.2/source/xdipo_rdc.f PARAxplor-nih-2.9.2/source/xdipo_rdc.f --- xplor-nih-2.9.2/source/xdipo_rdc.f 1970-01-01 01:00:00.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/xdipo_rdc.f 2003-12-05 14:34:36.000000000 +0100 @@ -0,0 +1,1216 @@ +C ------------------------------------------------------------------ + SUBROUTINE EXRDC (EDI,WHICH) +C ------------------------------------------------------------------ +C Calls EXRDC2, which does the actual energy calculation +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "xdipo_rdc.fcm" + INCLUDE "heap.fcm" + + DOUBLE PRECISION EDI + CHARACTER*7 WHICH + + CALL EXRDC2(EDI,HEAP(XRDCIPTR),HEAP(XRDCJPTR),HEAP(XRDCKPTR), + & HEAP(XRDCLPTR),HEAP(XRDCMPTR),HEAP(XRDCNPTR), + & HEAP(XRDCILST),HEAP(XRDCJLST),HEAP(XRDCKLST), + & HEAP(XRDCLLST),HEAP(XRDCMLST),HEAP(XRDCNLST), + & HEAP(XRDCOBSPTR),HEAP(XRDCERRPTR), + & HEAP(CALCXRDCPTR),WHICH) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE EXRDC2 (EDI,ATOMIPTR,ATOMJPTR,ATOMKPTR,ATOMLPTR, + & ATOMMPTR,ATOMNPTR,ATOMILST,ATOMJLST, + & ATOMKLST,ATOMLLST,ATOMMLST,ATOMNLST, + & XRDCOBS,XRDCERR,XRDCCALC,WHICH) +C ------------------------------------------------------------------ +C Calculates residual dipolar coupling energies +C +C Energies are of the form: +C E=K*(DELTARDC**2) +C where : +C K=FORCE CONSTANT +C DELTARDC=CALCULATED RDC-OBSERVED RDC +C and residual dipolar coupling function is defined as: +C RDC=A1*(3*COS(THETA)^2-1)+(3/2)*A2*(SIN(THETA)^2*COS(2*PHI)) +C +C WHICH is a flag that switches between energy/force calculations +C and RDC calculations for violations +C +C Note that: +C Atom I=Origin of the tensor +C Atom J=Z direction +C Atom K=X direction +C Atom L=Y direction +C Atom M=Coupled nucleus (e.g. N) +C Atom N=Detected nucleus (e.g. H) +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "psf.fcm" + INCLUDE "coord.fcm" + INCLUDE "numbers.fcm" + INCLUDE "deriv.fcm" + INCLUDE "xdipo_rdc.fcm" + INCLUDE "consta.fcm" + + DOUBLE PRECISION DIPOE + COMMON /DIPOENERGY/DIPOE + INTEGER ATOMIPTR(*),ATOMJPTR(*),ATOMKPTR(*) + INTEGER ATOMLPTR(*),ATOMMPTR(*),ATOMNPTR(*) + INTEGER ATOMILST(*),ATOMJLST(*),ATOMKLST(*) + INTEGER ATOMLLST(*),ATOMMLST(*),ATOMNLST(*) + DOUBLE PRECISION XRDCOBS(*),XRDCERR(*) + DOUBLE PRECISION XRDCCALC(*),EDI + CHARACTER*7 WHICH +C ------------------------------------------------------------------ +C Local variables +C ------------------------------------------------------------------ + INTEGER COUNT,CLASS,MM,NN,II,I,J,K,L,NA,M,N,SWAP + PARAMETER (NA=20) + INTEGER DMM,DNN,OUTFLAG(NA),MCOUNT,NCOUNT,PCOUNT + DOUBLE PRECISION XI,XJ,XK,XL,XM,XN,YI,YJ,YK,YL,YM,YN, + & ZI,ZJ,ZK,ZL,ZM,ZN, + & E,DF1(NA),DF2(NA),DF3(NA), + & PHI(NA),DPXI(NA),DPYI(NA),DPZI(NA), + & DPXJ(NA),DPYJ(NA),DPZJ(NA),DPXK(NA), + & DPYK(NA),DPZK(NA),DPXL(NA),DPYL(NA), + & DPZL(NA),DPXM(NA),DPYM(NA),DPZM(NA), + & DPXN(NA),DPYN(NA),DPZN(NA), + & THETA(NA),DTXI(NA),DTYI(NA),DTZI(NA), + & DTXJ(NA),DTYJ(NA),DTZJ(NA),DTXK(NA), + & DTYK(NA),DTZK(NA),DTXL(NA),DTYL(NA), + & DTZL(NA),DTXM(NA),DTYM(NA),DTZM(NA), + & DTXN(NA),DTYN(NA),DTZN(NA),DRXM(NA), + & DRYM(NA),DRZM(NA),DRXN(NA),DRYN(NA), + & DRZN(NA),DIST, + & O1,O2,O3,O4 ,O5 ,O6 ,O7 ,O8,O9,O10, + & O11,O12,O13,O14,O15,O16,O17,O18,O19, + & O20,O21,O22,O23,O24,O25,O26,O27,O28, + & O29,O30,O31,O32,O33,O34,O35,O36,O37, + & O38,O39,O40,O41,O42,O43,O44,O45,O46, + & O47,O48,O49,O50,O51,O52,O53,O54,O55, + & O56,O57,O58,O59,O60,O61,O62,O63,O64, + & O65,O66,O67,O68,O69,O70,O71,O72,O73, + & O74,O75,O76,O77,O78,O79,O80,O81,O82 + DOUBLE PRECISION OBSXRDC,ERRXRDC,K1, + & COEF1,COEF2,A,B,DELTA, + & DT,DP,DELTAXRDC, + & DD,DR,CALCXRDC +C ------------------------------------------------------------------ +C Zero out partial energy +C ------------------------------------------------------------------ + EDI=ZERO + + CLASS=1 + K1=XRDCFORCES(1) + COEF1=XRDCCOEF1(1) + COEF2=XRDCCOEF2(1) + + COUNT=0 + + DO WHILE (COUNT.LT.RDCNUM) + COUNT=COUNT+1 +C ------------------------------------------------------------------ +C Reset individual E to zero +C ------------------------------------------------------------------ + E=0 + + DO WHILE ((XRDCASSNDX(CLASS).LT.COUNT).OR. + & (XRDCASSNDX(CLASS).EQ.0)) + CLASS=CLASS+1 + END DO + + IF (XRDCASSNDX(CLASS).EQ.XRDCASSNDX(CLASS-1)) THEN + COUNT=COUNT-1 + ENDIF + + K1=XRDCFORCES(CLASS) + COEF1=XRDCCOEF1(CLASS) + COEF2=XRDCCOEF2(CLASS) +C ------------------------------------------------------------------ +C Note there should only be one atom for I, J, K, and L +C ------------------------------------------------------------------ + I=ATOMIPTR(COUNT)+1 + J=ATOMJPTR(COUNT)+1 + K=ATOMKPTR(COUNT)+1 + L=ATOMLPTR(COUNT)+1 + + XI=X(ATOMILST(I)) + XJ=X(ATOMJLST(J)) + XK=X(ATOMKLST(K)) + XL=X(ATOMLLST(L)) + YI=Y(ATOMILST(I)) + YJ=Y(ATOMJLST(J)) + YK=Y(ATOMKLST(K)) + YL=Y(ATOMLLST(L)) + ZI=Z(ATOMILST(I)) + ZJ=Z(ATOMJLST(J)) + ZK=Z(ATOMKLST(K)) + ZL=Z(ATOMLLST(L)) + + OBSXRDC=XRDCOBS(COUNT) + ERRXRDC=XRDCERR(COUNT) +C ------------------------------------------------------------------ +C Initialzie calculated residual dipolar coupling and counter +C ------------------------------------------------------------------ + CALCXRDC=0 + II=0 + MCOUNT=0 + NCOUNT=0 +C ------------------------------------------------------------------ +C Check for correct permutations of paired atoms to get the nucleus- +C nucleus vector. This depends solely on the order of the assigned +C atoms. OUTFLAG=1 indicates that the permutation is not allowed. +C ------------------------------------------------------------------ + DMM=ATOMMPTR(COUNT+1)-ATOMMPTR(COUNT) + DNN=ATOMNPTR(COUNT+1)-ATOMNPTR(COUNT) + DO MM=ATOMMPTR(COUNT)+1,ATOMMPTR(COUNT+1) + MCOUNT=MCOUNT+1 + NCOUNT=0 + DO NN=ATOMNPTR(COUNT)+1,ATOMNPTR(COUNT+1) + NCOUNT=NCOUNT+1 + II=II+1 + IF ((DMM.GT.1).AND.(DNN.GT.1)) THEN + IF (MCOUNT.EQ.NCOUNT) THEN + OUTFLAG(II)=0 + ELSE + OUTFLAG(II)=1 + ENDIF + ELSE + OUTFLAG(II)=0 + END IF + END DO + END DO + + II=0 + PCOUNT=0 + DO MM=ATOMMPTR(COUNT)+1,ATOMMPTR(COUNT+1) + DO NN=ATOMNPTR(COUNT)+1,ATOMNPTR(COUNT+1) + II=II+1 + IF (OUTFLAG(II).NE.1) THEN + PCOUNT=PCOUNT+1 + XM=X(ATOMMLST(MM)) + XN=X(ATOMNLST(NN)) + YM=Y(ATOMMLST(MM)) + YN=Y(ATOMNLST(NN)) + ZM=Z(ATOMMLST(MM)) + ZN=Z(ATOMNLST(NN)) +C ------------------------------------------------------------------ +C Calculate THETA, PHI and derivatives with respect to X, Y, Z of +C atoms I, J, K, L, M, N. This was done with Mathematica by Nico! +C ------------------------------------------------------------------ + O1=-XI+XJ + O2=O1**2 + O3=-YI+YJ + O4=O3**2 + O5=-ZI+ZJ + O6=O5**2 + O7=O2+O4+O6 + O8=SQRT(O7) + O9=1/O8 + O10=-XM+XN + O11=O1*O10 + O12=-YM+YN + O13=O12*O3 + O14=-ZM+ZN + O15=O14*O5 + O16=O11+O13+O15 + O17=O10**2 + O18=O12**2 + O19=O14**2 + O20=O17+O18+O19 + O21=SQRT(O20) + O22=1/O21 + O23=1/O7 + O24=O16**2 + O25=1/O20 + O26=O23*O24*O25 + O27=SQRT(1.D0-O26) + O28=1/O27 + O29=XM-XN + O30=O7**(3.D0/2.D0) + O31=1/O30 + O32=O1*O16*O22*O31 + O33=YM-YN + O34=O16*O22*O3*O31 + O35=ZM-ZN + O36=O16*O22*O31*O5 + O37=O20**(3.D0/2.D0) + O38=1/O37 + O39=O10*O16*O38*O9 + O40=O12*O16*O38*O9 + O41=O14*O16*O38*O9 + O42=-XI+XK + O43=O42**2 + O44=-YI+YK + O45=O44**2 + O46=-ZI+ZK + O47=O46**2 + O48=O43+O45+O47 + O49=SQRT(O48) + O50=-XI+XL + O51=O50**2 + O52=-YI+YL + O53=O52**2 + O54=-ZI+ZL + O55=O54**2 + O56=O51+O53+O55 + O57=SQRT(O56) + O58=1/O57 + O59=O10*O42 + O60=O12*O44 + O61=O14*O46 + O62=O59+O60+O61 + O63=1/O62 + O64=O10*O50 + O65=O12*O52 + O66=O14*O54 + O67=O64+O65+O66 + O68=O62**2 + O69=1/O68 + O70=O56**(3.D0/2.D0) + O71=1/O70 + O72=O49*O50*O63*O67*O71 + O73=1/O49 + O74=O42*O58*O63*O67*O73 + O75=1/O56 + O76=O67**2 + O77=O48*O69*O75*O76 + O78=1/(1.D0+O77) + O79=O49*O52*O63*O67*O71 + O80=O44*O58*O63*O67*O73 + O81=O49*O54*O63*O67*O71 + O82=O46*O58*O63*O67*O73 + + THETA(II)=ACOS(O16*O22*O9) + + DTXI(II)=-(O28*(O32+O22*O29*O9)) + DTYI(II)=-(O28*(O34+O22*O33*O9)) + DTZI(II)=-(O28*(O36+O22*O35*O9)) + DTXJ(II)=-(O28*(-O32+O10*O22*O9)) + DTYJ(II)=-(O28*(-O34+O12*O22*O9)) + DTZJ(II)=-(O28*(-O36+O14*O22*O9)) + DTXK(II)=0 + DTYK(II)=0 + DTZK(II)=0 + DTXL(II)=0 + DTYL(II)=0 + DTZL(II)=0 + DTXM(II)=-(O28*(O39+O22*O9*(XI-XJ))) + DTYM(II)=-(O28*(O40+O22*O9*(YI-YJ))) + DTZM(II)=-(O28*(O41+O22*O9*(ZI-ZJ))) + DTXN(II)=-(O28*(-O39+O1*O22*O9)) + DTYN(II)= -(O28*(-O40+O22*O3*O9)) + DTZN(II)=-(O28*(-O41+O22*O5*O9)) + + PHI(II)=ATAN(O49*O58*O63*O67) + + DPXI(II)=O78*(O29*O49*O58*O63-O29*O49*O58*O67*O69+ + & O72-O74) + DPYI(II)=O78*(O33*O49*O58*O63-O33*O49*O58*O67*O69+ + & O79-O80) + DPZI(II)=O78*(O35*O49*O58*O63-O35*O49*O58*O67*O69+ + & O81-O82) + DPXJ(II)=0 + DPYJ(II)=0 + DPZJ(II)=0 + DPXK(II)=(-(O10*O49*O58*O67*O69)+O74)*O78 + DPYK(II)=O78*(-(O12*O49*O58*O67*O69)+O80) + DPZK(II)=O78*(-(O14*O49*O58*O67*O69)+O82) + DPXL(II)=(O10*O49*O58*O63-O72)*O78 + DPYL(II)=O78*(O12*O49*O58*O63-O79) + DPZL(II)=O78*(O14*O49*O58*O63-O81) + DPXM(II)=O78*(-(O49*O58*O67*O69*(XI-XK))+ + & O49*O58*O63*(XI-XL)) + DPYM(II)=O78*(-(O49*O58*O67*O69*(YI-YK))+ + & O49*O58*O63*(YI-YL)) + DPZM(II)=O78*(-(O49*O58*O67*O69*(ZI-ZK))+ + & O49*O58*O63*(ZI-ZL)) + DPXN(II)=(O49*O50*O58*O63-O42*O49*O58*O67*O69)*O78 + DPYN(II)=(O49*O52*O58*O63-O44*O49*O58*O67*O69)*O78 + DPZN(II)=(O49*O54*O58*O63-O46*O49*O58*O67*O69)*O78 +C ------------------------------------------------------------------ +C Calculate the distance between the two atoms M, N (coupled nuclei) +C and calculate the derivatives as well +C ------------------------------------------------------------------ + O1=-XM+XN + O2=O1**2 + O3=-YM+YN + O4=O3**2 + O5=-ZM+ZN + O6=O5**2 + O7=O2+O4+O6 + O8=SQRT(O7) + O9=1/O8 + O10=O1*O9 + O11=O3*O9 + O12=O5*O9 + DIST=O8 + DRXM(II)=-O10 + DRYM(II)=-O11 + DRZM(II)=-O12 + DRXN(II)=O10 + DRYN(II)=O11 + DRZN(II)=O12 +C ------------------------------------------------------------------ +C Calculate the residual dipolar coupling +C ------------------------------------------------------------------ + A=COEF1 + B=COEF2 + CALCXRDC=A*(3.*COS(THETA(II))*COS(THETA(II))-1.)+ + & B*(3./2.)*(SIN(THETA(II))*SIN(THETA(II))* + & COS(2.*PHI(II))) + END IF + END DO + END DO + + IF (WHICH.EQ.'ANALYZE') THEN + XRDCCALC(COUNT)=CALCXRDC + END IF + + DELTA=CALCXRDC-OBSXRDC +C ------------------------------------------------------------------ +C Adjust the deviation based on the error range +C ------------------------------------------------------------------ + IF ((DELTA.LT.0.000).AND.(ABS(DELTA).GT.ERRXRDC)) THEN + DELTAXRDC=DELTA+ERRXRDC + ELSE IF ((DELTA.GT.0.000).AND.(DELTA.GT.ERRXRDC)) THEN + DELTAXRDC=DELTA-ERRXRDC + ELSE + DELTAXRDC=0.0 + END IF + + DD=DELTAXRDC + + II=0 + PCOUNT=0 + DO MM=ATOMMPTR(COUNT)+1,ATOMMPTR(COUNT+1) + DO NN=ATOMNPTR(COUNT)+1,ATOMNPTR(COUNT+1) + II=II+1 + IF (OUTFLAG(II).NE.1) THEN + PCOUNT=PCOUNT+1 +C ------------------------------------------------------------------ +C Taking care of derivative +C ------------------------------------------------------------------ + DT=(A*(-6.)*SIN(THETA(II))* + & COS(THETA(II))+B*3.*SIN(THETA(II))* + & COS(THETA(II))*COS(2.*PHI(II))) + DP=(B*(-3.)*SIN(THETA(II))* + & SIN(THETA(II))*SIN(2.*PHI(II))) + DR=0.0 + IF (DD.EQ.0.0) THEN + E=E+0.0 + DF1(II)=0.0 + DF2(II)=0.0 + DF3(II)=0.0 + ELSE + E=E+K1*(DELTAXRDC**2) + DF1(II)=2*K1*DELTAXRDC*DT + DF2(II)=2*K1*DELTAXRDC*DP + DF3(II)=2*K1*DELTAXRDC*DR + END IF + END IF + END DO + END DO +C ------------------------------------------------------------------ +C Accumulate energy +C ------------------------------------------------------------------ + EDI=EDI+E +C ------------------------------------------------------------------ +C Now update forces if in energy/force mode +C ------------------------------------------------------------------ + IF (WHICH.NE.'ANALYZE') THEN + II=0 + DO MM=ATOMMPTR(COUNT)+1,ATOMMPTR(COUNT+1) + DO NN=ATOMNPTR(COUNT)+1,ATOMNPTR(COUNT+1) + II=II+1 + IF (OUTFLAG(II).NE.1) THEN + DX(ATOMILST(I))=DX(ATOMILST(I))+ + & DF1(II)*DTXI(II)+DF2(II)*DPXI(II) + DY(ATOMILST(I))=DY(ATOMILST(I))+ + & DF1(II)*DTYI(II)+DF2(II)*DPYI(II) + DZ(ATOMILST(I))=DZ(ATOMILST(I))+ + & DF1(II)*DTZI(II)+DF2(II)*DPZI(II) + DX(ATOMJLST(J))=DX(ATOMJLST(J))+ + & DF1(II)*DTXJ(II)+DF2(II)*DPXJ(II) + DY(ATOMJLST(J))=DY(ATOMJLST(J))+ + & DF1(II)*DTYJ(II)+DF2(II)*DPYJ(II) + DZ(ATOMJLST(J))=DZ(ATOMJLST(J))+ + & DF1(II)*DTZJ(II)+DF2(II)*DPZJ(II) + DX(ATOMKLST(K))=DX(ATOMKLST(K))+ + & DF1(II)*DTXK(II)+DF2(II)*DPXK(II) + DY(ATOMKLST(K))=DY(ATOMKLST(K))+ + & DF1(II)*DTYK(II)+DF2(II)*DPYK(II) + DZ(ATOMKLST(K))=DZ(ATOMKLST(K))+ + & DF1(II)*DTZK(II)+DF2(II)*DPZK(II) + DX(ATOMLLST(L))=DX(ATOMLLST(L))+ + & DF1(II)*DTXL(II)+DF2(II)*DPXL(II) + DY(ATOMLLST(L))=DY(ATOMLLST(L))+ + & DF1(II)*DTYL(II)+DF2(II)*DPYL(II) + DZ(ATOMLLST(L))=DZ(ATOMLLST(L))+ + & DF1(II)*DTZL(II)+DF2(II)*DPZL(II) + DX(ATOMMLST(MM))=DX(ATOMMLST(MM))+ + & DF1(II)*DTXM(II)+DF2(II)* + & DPXM(II)+DF3(II)*DRXM(II) + DY(ATOMMLST(MM))=DY(ATOMMLST(MM))+ + & DF1(II)*DTYM(II)+DF2(II)* + & DPYM(II)+DF3(II)*DRYM(II) + DZ(ATOMMLST(MM))=DZ(ATOMMLST(MM))+ + & DF1(II)*DTZM(II)+DF2(II)* + & DPZM(II)+DF3(II)*DRZM(II) + DX(ATOMNLST(NN))=DX(ATOMNLST(NN))+ + & DF1(II)*DTXN(II)+DF2(II)* + & DPXN(II)+DF3(II)*DRXN(II) + DY(ATOMNLST(NN))=DY(ATOMNLST(NN))+ + & DF1(II)*DTYN(II)+DF2(II)* + & DPYN(II)+DF3(II)*DRYN(II) + DZ(ATOMNLST(NN))=DZ(ATOMNLST(NN))+ + & DF1(II)*DTZN(II)+DF2(II)* + & DPZN(II)+DF3(II)*DRZN(II) + END IF + END DO + END DO + END IF + END DO + + DIPOE=EDI + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE READXRDC +C ------------------------------------------------------------------ +C Reads in residual dipolar coupling information +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "comand.fcm" + INCLUDE "xdipo_rdc.fcm" + INCLUDE "funct.fcm" + INCLUDE "psf.fcm" + INCLUDE "heap.fcm" + INCLUDE "numbers.fcm" + + LOGICAL FDSWITCH,FDSAVE,FDERR + COMMON /FDPARAM/FDSWITCH,FDSAVE,FDERR + DOUBLE PRECISION FDEPERC + INTEGER FDENUMC + COMMON /FDERRORE/FDEPERC,FDENUMC + + INTEGER COUNT,SPTR,OLDCLASS,OLDMAXXRDC,FNCLASS,FNUMC + DOUBLE PRECISION K1,K2,CUTOFF,COEF1,COEF2 + DOUBLE PRECISION FMEDPERC,FPERC + CHARACTER*4 THENAME + CHARACTER*132 NOMEFILE + INTEGER TOLSWI + + NOMEFILE='XDIPO_RDC.OUT' + + SPTR=ALLHP(INTEG4(NATOM)) + CALL PUSEND('XDIPO_RDC>') +862 CONTINUE + CALL NEXTWD('XDIPO_RDC>') + CALL MISCOM('XDIPO_RDC>',USED) + IF (.NOT.USED) THEN +C ------------------------------------------------------------------ +C Documentation +C ------------------------------------------------------------------ + IF (WD(1:4).EQ.'HELP') THEN + WRITE(DUNIT,'(10X,A)') + & ' XRDC {} END ', + & ' :== ', + & ' ASSIgn ', + & ' ', + & ' * Restraint: Metal Z X Y Coupled_nucleus', + & ' Detected_nucleus RDC Err *', + & ' CLASsification ', + & ' * Starts a new class *', + & ' TOLL <1|0>', + & ' * Switch tolerance in FANTALIN (1=Yes, 0=No)', + & ' COEFficient ', + & ' * Tensor parameters: A1 A2 *', + & ' FORCeconstant ', + & ' * Force constant for the current class *', + & ' NREStraints ', + & ' * Number of slots for restraints to allocate *', + & ' PRINt THREshold ', + & ' * Prints violations larger than the THREshold *', + & ' RESEt', + & ' * Erases the restraint table, but keeps NRES *', + & ' SAVE', + & ' * Filename to save RDC values *', + & ' FMED', + & ' * Calculates average tensor parameters *', + & ' ERRON', + & ' * Switches Monte Carlo error evaluation ON *', + & ' ERROFF', + & ' * Switches Monte Carlo error evaluation OFF *', + & ' FON', + & ' * Switches tensor auto-update mode ON *', + & ' FOFF', + & ' * Switches tensor auto-update mode OFF *', + & ' SON', + & ' * Switches saving mode ON *', + & ' SOFF', + & ' * Switches saving mode OFF *', + & ' FRUN', + & ' * Runs FANTALIN *' +C ------------------------------------------------------------------ +C About FANTALIN +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'TOLL') THEN + CALL NEXTI('Tolerance setting in FANTALIN (1/0) =',TOLSWI) + IF (TOLSWI.EQ.1) THEN + DIPFLG=.TRUE. + WRITE (DUNIT,'(A)') 'Tolerance ON in FANTALIN.' + ELSE IF (TOLSWI.EQ.0) THEN + DIPFLG=.FALSE. + WRITE (DUNIT,'(A)') 'Tolerance OFF in FANTALIN.' + ELSE + WRITE (DUNIT,'(A)') 'Unknown switch. Using default.' + END IF + ELSE IF (WD(1:4).EQ.'SAVE') THEN + CALL NEXTFI('Filename =',NOMEFILE) + PRINT*,'Name of the file to save RDC values :',NOMEFILE + ELSE IF (WD(1:3).EQ.'FME') THEN + CALL NEXTF('Percentage for average (0/100) =',FMEDPERC) + CALL NEXTI('Average on class number =',FNCLASS) + IF (FNCLASS.GT.NCLASSES) THEN + PRINT*,'%FRUN-ERR: This class does not exist...' + ELSE + CALL FANTADIPOMEDIA(FNCLASS,FMEDPERC) + END IF + ELSE IF (WD(1:5).EQ.'ERRON') THEN + CALL NEXTI('Number of cycles =',FDENUMC) + CALL NEXTF('Percentage to discard (0/100) =',FDEPERC) + FDEPERC=FDEPERC/100.0 + FDERR=1 + PRINT*,'Monte Carlo error evaluation ON' + ELSE IF (WD(1:6).EQ.'ERROFF') THEN + FDERR=0 + PRINT*,'Monte Carlo error evaluation OFF' + ELSE IF (WD(1:3).EQ.'FON') THEN + FDSWITCH=1 + PRINT*,'Tensor auto-update mode ON' + ELSE IF (WD(1:4).EQ.'FOFF') THEN + FDSWITCH=0 + PRINT*,'Tensor auto-update mode OFF' + ELSE IF (WD(1:3).EQ.'SON') THEN + FDSAVE=1 + PRINT*,'Saving mode ON' + ELSE IF (WD(1:4).EQ.'SOFF') THEN + FDSAVE=0 + PRINT*,'Saving mode OFF' + ELSE IF (WD(1:4).EQ.'FRUN') THEN + CALL NEXTI('FANTALIN on class number =',FNCLASS) + IF (FNCLASS.GT.NCLASSES.OR.FNCLASS.EQ.0) THEN + PRINT*,'%FRUN-ERR: This class does not exist...' + ELSE + CALL FANTADIPO(HEAP(XRDCIPTR),HEAP(XRDCJPTR), + & HEAP(XRDCKPTR),HEAP(XRDCLPTR), + & HEAP(XRDCMPTR),HEAP(XRDCNPTR), + & HEAP(XRDCILST),HEAP(XRDCJLST), + & HEAP(XRDCKLST),HEAP(XRDCLLST), + & HEAP(XRDCMLST),HEAP(XRDCNLST), + & HEAP(XRDCOBSPTR),HEAP(XRDCERRPTR), + & FNCLASS,NOMEFILE) + END IF +C ------------------------------------------------------------------ +C Get class name. Determine if it's an already-defined class. +C Insert a new class if it's not. +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'CLAS') THEN + OLDCLASS=CURCLASS + CALL NEXTA4('Class name =',THENAME) + MODE=NEW + DO COUNT=1,NCLASSES + IF (XRDCCLASSNAMES(COUNT).EQ.THENAME) THEN + MODE=UPDATE + CURCLASS=COUNT + END IF + END DO + IF (MODE.EQ.NEW) THEN +C ------------------------------------------------------------------ +C Make sure you can't add more than the maximum number of classes +C ------------------------------------------------------------------ + IF (OLDCLASS.EQ.MAXXRDCCLASSES) THEN + CALL DSPERR('XDIPO_RDC','Too many classes.') + CALL DSPERR('XDIPO_RDC', + & 'Increase MAXXRDCCLASSES and recompile.') + CALL WRNDIE(-5, 'READXRDC', + & 'Too many anisotropy classes.') + END IF + NCLASSES=NCLASSES+1 + CURCLASS=NCLASSES + XRDCCLASSNAMES(CURCLASS)=THENAME +C ------------------------------------------------------------------ +C If this isn't the first class, close off the old class +C ------------------------------------------------------------------ + IF (NCLASSES.GT.1) THEN + XRDCASSNDX(OLDCLASS)=RDCNUM + END IF + END IF +C ------------------------------------------------------------------ +C Set force constant for current class +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'FORC') THEN +C ------------------------------------------------------------------ +C Start a default class if there isn't one defined +C ------------------------------------------------------------------ + IF (CURCLASS.EQ.0) THEN + NCLASSES=1 + CURCLASS=1 + END IF + CALL NEXTF('Force constant =',K1) + WRITE(DUNIT,'(A,A,A,F8.3)') + & 'Setting force constant for class ', + & XRDCCLASSNAMES(CURCLASS),' to ',K1 + XRDCFORCES(CURCLASS)=K1 +C ------------------------------------------------------------------ +C Set coefficient constant for current class +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'COEF') THEN + CALL NEXTF('Coefficient A1 =',COEF1) + CALL NEXTF('Coefficient A2 =',COEF2) +C ------------------------------------------------------------------ +C Start a default class if there isn't one defined +C ------------------------------------------------------------------ + IF (CURCLASS.EQ.0) THEN + NCLASSES=1 + CURCLASS=1 + END IF + WRITE(DUNIT,'(A,A,A,F8.3,F8.3)') + & 'Setting coefficients for class ', + & XRDCCLASSNAMES(CURCLASS),' to ',COEF1,COEF2 + XRDCCOEF1(CURCLASS)=COEF1 + XRDCCOEF2(CURCLASS)=COEF2 +C ------------------------------------------------------------------ +C Reset residual dipolar couplings database +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'RESE') THEN + CALL XRDCDEFAULTS + CALL ALLOCXRDC(0,MAXXRDC) +C ------------------------------------------------------------------ +C Change number of assignment slots +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'NRES') THEN + OLDMAXXRDC=MAXXRDC + CALL NEXTI('Number of slots =',MAXXRDC) + CALL ALLOCXRDC(OLDMAXXRDC,MAXXRDC) + WRITE(DUNIT,'(A,I8,A)') + & 'XDIPO_RDC: Allocating space for',MAXXRDC, + & 'number of restraints.' +C ------------------------------------------------------------------ +C Read in an assignment +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'ASSI') THEN +C ------------------------------------------------------------------ +C Make sure you can't add more assignments than you have slots for +C ------------------------------------------------------------------ + IF (XRDCMX.EQ.MAXXRDC) THEN + CALL DSPERR('XDIPO_RDC','Too many assignments.') + CALL DSPERR('XDIPO_RDC', + & 'Increasing NREStraints by 100.') + OLDMAXXRDC=MAXXRDC + MAXXRDC=MAXXRDC+100 + CALL ALLOCXRDC(OLDMAXXRDC,MAXXRDC) + END IF +C ------------------------------------------------------------------ +C If there isn't a class specified, start a default class +C ------------------------------------------------------------------ + IF (CURCLASS.EQ.0) THEN + NCLASSES=1 + CURCLASS=1 + END IF + CALL READXRDC2(HEAP(XRDCIPTR),HEAP(XRDCJPTR), + & HEAP(XRDCKPTR),HEAP(XRDCLPTR), + & HEAP(XRDCMPTR),HEAP(XRDCNPTR), + & HEAP(XRDCILST),HEAP(XRDCJLST), + & HEAP(XRDCKLST),HEAP(XRDCLLST), + & HEAP(XRDCMLST),HEAP(XRDCNLST), + & HEAP(XRDCOBSPTR),HEAP(XRDCERRPTR), + & HEAP(SPTR)) +C ------------------------------------------------------------------ +C Print violations +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'PRIN') THEN + CALL NEXTWD('PRINT>') + IF (WD(1:4).NE.'THRE') THEN + CALL DSPERR('XDIPO_RDC', + & 'PRINt expects THREshold parameter.') + ELSE + CALL NEXTF('THREshold =',CUTOFF) + IF (CUTOFF.LT.ZERO) THEN + CALL DSPERR('XDIPO_RDC', + & 'Cutoff must be positive.') + CUTOFF=ABS(CUTOFF) + END IF + CALL NEXTA4('ALL OR CLASs>',THENAME) + IF (THENAME(1:3).EQ.'ALL') THEN + DO COUNT=1,NCLASSES + PRINTCLASS(COUNT)=.TRUE. + END DO + ELSE IF (THENAME(1:4).EQ.'CLAS') THEN + CALL NEXTA4('CLASS NAME =',THENAME) + DO COUNT=1,NCLASSES + IF (XRDCCLASSNAMES(COUNT).EQ.THENAME) THEN + PRINTCLASS(COUNT)=.TRUE. + ELSE + PRINTCLASS(COUNT)=.FALSE. + END IF + END DO + ELSE + DO COUNT=1,NCLASSES + PRINTCLASS(COUNT)=.TRUE. + END DO + END IF + CALL PRINTXRDC(CUTOFF,HEAP(CALCXRDCPTR), + & HEAP(XRDCOBSPTR),HEAP(XRDCERRPTR), + & HEAP(XRDCIPTR),HEAP(XRDCJPTR), + & HEAP(XRDCKPTR),HEAP(XRDCLPTR), + & HEAP(XRDCMPTR),HEAP(XRDCNPTR), + & HEAP(XRDCILST),HEAP(XRDCJLST), + & HEAP(XRDCKLST),HEAP(XRDCLLST), + & HEAP(XRDCMLST),HEAP(XRDCNLST)) + END IF +C ------------------------------------------------------------------ +C Check for END statement +C ------------------------------------------------------------------ + ELSE + CALL CHKEND('XDIPO_RDC>',DONE) + END IF + END IF + IF (.NOT.DONE) GOTO 862 + DONE=.FALSE. + CALL FREHP(SPTR,INTEG4(NATOM)) + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE ALLOCXRDC (OLDSIZE,NEWSIZE) +C ------------------------------------------------------------------ +C Reset residual dipolar coupling arrays to hold size entries + +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + + IMPLICIT NONE + + INCLUDE "funct.fcm" + INCLUDE "xdipo_rdc.fcm" + + INTEGER OLDSIZE, NEWSIZE + + IF (OLDSIZE.NE.0) THEN + CALL FREHP(XRDCIPTR,INTEG4(OLDSIZE)) + CALL FREHP(XRDCJPTR,INTEG4(OLDSIZE)) + CALL FREHP(XRDCKPTR,INTEG4(OLDSIZE)) + CALL FREHP(XRDCLPTR,INTEG4(OLDSIZE)) + CALL FREHP(XRDCMPTR,INTEG4(OLDSIZE)) + CALL FREHP(XRDCNPTR,INTEG4(OLDSIZE)) + + CALL FREHP(XRDCILST,INTEG4(OLDSIZE)) + CALL FREHP(XRDCJLST,INTEG4(OLDSIZE)) + CALL FREHP(XRDCKLST,INTEG4(OLDSIZE)) + CALL FREHP(XRDCLLST,INTEG4(OLDSIZE)) + CALL FREHP(XRDCMLST,INTEG4(OLDSIZE)) + CALL FREHP(XRDCNLST,INTEG4(OLDSIZE)) + + CALL FREHP(XRDCOBSPTR,IREAL8(OLDSIZE)) + CALL FREHP(XRDCERRPTR,IREAL8(OLDSIZE)) + CALL FREHP(CALCXRDCPTR,IREAL8(OLDSIZE)) + END IF + XRDCIPTR=ALLHP(INTEG4(NEWSIZE)) + XRDCJPTR=ALLHP(INTEG4(NEWSIZE)) + XRDCKPTR=ALLHP(INTEG4(NEWSIZE)) + XRDCLPTR=ALLHP(INTEG4(NEWSIZE)) + XRDCMPTR=ALLHP(INTEG4(NEWSIZE)) + XRDCNPTR=ALLHP(INTEG4(NEWSIZE)) + + XRDCILST=ALLHP(INTEG4(NEWSIZE)) + XRDCJLST=ALLHP(INTEG4(NEWSIZE)) + XRDCKLST=ALLHP(INTEG4(NEWSIZE)) + XRDCLLST=ALLHP(INTEG4(NEWSIZE)) + XRDCMLST=ALLHP(INTEG4(NEWSIZE)) + XRDCNLST=ALLHP(INTEG4(NEWSIZE)) + + XRDCOBSPTR=ALLHP(IREAL8(NEWSIZE)) + XRDCERRPTR=ALLHP(IREAL8(NEWSIZE)) + CALCXRDCPTR=ALLHP(IREAL8(NEWSIZE)) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE XRDCDEFAULTS +C ------------------------------------------------------------------ +C Sets up defaults +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "xdipo_rdc.fcm" + + INTEGER COUNT + + MODE=NEW + MAXXRDC=200 + XRDCMX=200 + RDCNUM=0 + NCLASSES=0 + CURCLASS=0 + DIPFLG=.FALSE. + DO COUNT=1, MAXXRDCCLASSES + XRDCCLASSNAMES(COUNT)='DEFAULT' + XRDCASSNDX(COUNT)=0 + XRDCFORCES(COUNT)=5.0 + XRDCCOEF1(COUNT)=1.0 + XRDCCOEF2(COUNT)=1.0 + END DO + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE READXRDC2(ATOMIPTR,ATOMJPTR,ATOMKPTR,ATOMLPTR,ATOMMPTR, + & ATOMNPTR,ATOMILST,ATOMJLST,ATOMKLST,ATOMLLST, + & ATOMMLST,ATOMNLST,XRDCOBS,XRDCERR,SEL) +C ------------------------------------------------------------------ +C Reads actual residual dipolar coupling assignments into arrays +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "psf.fcm" + INCLUDE "coord.fcm" + INCLUDE "xdipo_rdc.fcm" + + INTEGER ATOMIPTR(*),ATOMJPTR(*),ATOMKPTR(*) + INTEGER ATOMLPTR(*),ATOMMPTR(*),ATOMNPTR(*) + INTEGER ATOMILST(*),ATOMJLST(*),ATOMKLST(*) + INTEGER ATOMLLST(*),ATOMMLST(*),ATOMNLST(*),SEL(*) + INTEGER ITMP,JTMP,KTMP,LTMP,II + DOUBLE PRECISION XRDCOBS(*),XRDCERR(*) +C ------------------------------------------------------------------ +C Local variables +C ------------------------------------------------------------------ + INTEGER NSEL,INSERTPOS,COUNT,CURSTOP,OTHERSTOP,NFLAG + INTEGER I + DOUBLE PRECISION XRDCO,XRDCE +C ------------------------------------------------------------------ +C If we're in UPDATE mode, make a space for the new line +C ------------------------------------------------------------------ + NFLAG=0 + IF (MODE.EQ.UPDATE) THEN + DO COUNT=RDCNUM+1,XRDCASSNDX(CURCLASS)+1, -1 + ATOMIPTR(COUNT)=ATOMIPTR(COUNT-1) + ATOMJPTR(COUNT)=ATOMJPTR(COUNT-1) + ATOMKPTR(COUNT)=ATOMKPTR(COUNT-1) + ATOMLPTR(COUNT)=ATOMLPTR(COUNT-1) + ATOMMPTR(COUNT)=ATOMMPTR(COUNT-1) + ATOMNPTR(COUNT)=ATOMNPTR(COUNT-1) + XRDCOBS(COUNT)=XRDCOBS(COUNT-1) + XRDCERR(COUNT)=XRDCERR(COUNT-1) + END DO + CURSTOP=XRDCASSNDX(CURCLASS) + DO COUNT=1,NCLASSES + OTHERSTOP=XRDCASSNDX(COUNT) + IF (OTHERSTOP.GT.CURSTOP) THEN + XRDCASSNDX(COUNT)=OTHERSTOP+1 + END IF + END DO + XRDCASSNDX(CURCLASS)=CURSTOP+1 + INSERTPOS=CURSTOP + RDCNUM=RDCNUM+1 + ELSE + RDCNUM=RDCNUM+1 + INSERTPOS=RDCNUM + XRDCASSNDX(CURCLASS)=INSERTPOS + END IF +C ------------------------------------------------------------------ +C Reading in the atom selection in the restraint table +C ------------------------------------------------------------------ + CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.) + IF (NSEL.GT.1) THEN + CALL DSPERR('XDIPO_RDC', + & 'More than 1 atom in SEL for atom I. Using first') + END IF + IF (NSEL.EQ.0) THEN + CALL DSPERR('XDIPO_RDC','Error with atom selection') + NFLAG=1 + END IF + CALL MAKIND(SEL,NATOM,NSEL) + IF (INSERTPOS.EQ.1) ATOMIPTR(INSERTPOS)=0 + II=ATOMIPTR(INSERTPOS) + II=II+1 + IF (II.GT.XRDCMX) XRDCMX=II + ATOMILST(II)=SEL(1) + ATOMIPTR(INSERTPOS+1)=II + + CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.) + IF (NSEL.GT.1) THEN + CALL DSPERR('XDIPO_RDC', + & 'More than 1 atom in SEL for atom J. Using first') + END IF + IF (NSEL.EQ.0) THEN + CALL DSPERR('XDIPO_RDC','Error with atom selection') + NFLAG=1 + END IF + CALL MAKIND(SEL,NATOM,NSEL) + IF (INSERTPOS.EQ.1) ATOMJPTR(INSERTPOS)=0 + II=ATOMJPTR(INSERTPOS) + II=II+1 + IF (II.GT.XRDCMX) XRDCMX=II + ATOMJLST(II)=SEL(1) + ATOMJPTR(INSERTPOS+1)=II + + CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.) + IF (NSEL.GT.1) THEN + CALL DSPERR('XDIPO_RDC', + & 'More than 1 atom in SEL for atom K. Using first') + END IF + IF (NSEL.EQ.0) THEN + CALL DSPERR('XDIPO_RDC','Error with atom selection') + NFLAG=1 + END IF + CALL MAKIND(SEL, NATOM, NSEL) + IF (INSERTPOS.EQ.1) ATOMKPTR(INSERTPOS)=0 + II=ATOMKPTR(INSERTPOS) + II=II+1 + IF (II.GT.XRDCMX) XRDCMX=II + ATOMKLST(II)=SEL(1) + ATOMKPTR(INSERTPOS+1)=II + + CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.) + IF (NSEL.GT.1) THEN + CALL DSPERR('XDIPO_RDC', + & 'More than 1 atom in SEL for atom L. Using first') + END IF + IF (NSEL.EQ.0) THEN + CALL DSPERR('XDIPO_RDC','Error with atom selection') + NFLAG=1 + END IF + CALL MAKIND(SEL,NATOM,NSEL) + IF (INSERTPOS.EQ.1) ATOMLPTR(INSERTPOS)=0 + II=ATOMLPTR(INSERTPOS) + II=II+1 + IF (II.GT.XRDCMX) XRDCMX=II + ATOMLLST(II)=SEL(1) + ATOMLPTR(INSERTPOS+1)=II +C ------------------------------------------------------------------ +C Next two selections define the vector of the coupled nuclei. Note +C that the order of atoms picked is important. +C ------------------------------------------------------------------ + CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.) + IF (NSEL.GT.1) THEN + CALL DSPERR('XDIPO_RDC', + & 'More than 1 atom in SEL for atom M. Using first') + END IF + IF (NSEL.EQ.0) THEN + CALL DSPERR('XDIPO_RDC','Error with atom selection') + NFLAG=1 + END IF + CALL MAKIND(SEL,NATOM,NSEL) + IF (INSERTPOS.EQ.1) ATOMMPTR(INSERTPOS)=0 + II=ATOMMPTR(INSERTPOS) + II=II+1 + IF (II.GT.XRDCMX) XRDCMX=II + ATOMMLST(II)=SEL(1) + ATOMMPTR(INSERTPOS+1)=II + + CALL SELCTA(SEL,NSEL,X,Y,Z,.TRUE.) + IF (NSEL.GT.1) THEN + CALL DSPERR('XDIPO_RDC', + & 'More than 1 atom in SEL for atom N. Using first') + END IF + IF (NSEL.EQ.0) THEN + CALL DSPERR('XDIPO_RDC','Error with atom selection') + NFLAG=1 + END IF + CALL MAKIND(SEL,NATOM,NSEL) + IF (INSERTPOS.EQ.1) ATOMNPTR(INSERTPOS)=0 + II=ATOMNPTR(INSERTPOS) + II=II+1 + IF (II.GT.XRDCMX) XRDCMX=II + ATOMNLST(II)=SEL(1) + ATOMNPTR(INSERTPOS+1)=II +C ------------------------------------------------------------------ +C Reading in the observed residual dipolar coupling +C ------------------------------------------------------------------ + CALL NEXTF('Observed RDC =',XRDCO) + CALL NEXTF('Error in RDC =',XRDCE) + + XRDCOBS(INSERTPOS)=XRDCO + XRDCERR(INSERTPOS)=XRDCE +C ------------------------------------------------------------------ +C Check for error atom selection. If there is one then reset the +C counter for restraint +C ------------------------------------------------------------------ + IF (NFLAG.EQ.1) THEN + RDCNUM=RDCNUM-1 + END IF + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE XRDCINIT +C ------------------------------------------------------------------ +C Initializes residual dipolar coupling stuff +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "xdipo_rdc.fcm" + + CALL XRDCDEFAULTS + CALL ALLOCXRDC(0,MAXXRDC) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE PRINTXRDC (CUTOFF,XRDCCALC,XRDCOBS,XRDCERR, + & ATOMIPTR,ATOMJPTR,ATOMKPTR, + & ATOMLPTR,ATOMMPTR,ATOMNPTR, + & ATOMILST,ATOMJLST,ATOMKLST, + & ATOMLLST,ATOMMLST,ATOMNLST) +C ------------------------------------------------------------------ +C Prints residual dipolar couplings with deviation larger than cutoff, +C calculates RMSD and puts it into $RESULT +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "xdipo_rdc.fcm" + INCLUDE "comand.fcm" + INCLUDE "psf.fcm" + INCLUDE "numbers.fcm" + + DOUBLE PRECISION CUTOFF,XRDCCALC(*),XRDCOBS(*),XRDCERR(*) + INTEGER ATOMILST(*),ATOMJLST(*),ATOMKLST(*) + INTEGER ATOMLLST(*),ATOMMLST(*),ATOMNLST(*) + INTEGER ATOMIPTR(*),ATOMJPTR(*),ATOMKPTR(*) + INTEGER ATOMLPTR(*),ATOMMPTR(*),ATOMNPTR(*) +C ------------------------------------------------------------------ +C Local variables +C ------------------------------------------------------------------ + DOUBLE PRECISION CALCXRDC,OBSXRDC,DELTAXRDC,DELTA,DP + INTEGER COUNT,CLASS,I,J,K,L,M,N,NUM,II + DOUBLE PRECISION RMS,VIOLS,ERRXRDC,XRDCENERGY + DOUBLE COMPLEX DUMMY2 + LOGICAL PRINTTHISCLASS + + RMS=ZERO + VIOLS=ZERO + NUM=0 +C ------------------------------------------------------------------ +C Make sure that the array of calculated RDC is up to date +C ------------------------------------------------------------------ + CALL EXRDC(XRDCENERGY,'ANALYZE') + WRITE (PUNIT,'(A)') 'The following RDC''s have' + WRITE (PUNIT,'(A)') 'deviation larger than or' + WRITE (PUNIT,'(A)') 'equal to the cutoff:' +C ------------------------------------------------------------------ +C Write out first class heading +C ------------------------------------------------------------------ + CLASS=1 + PRINTTHISCLASS=PRINTCLASS(CLASS) + IF (PRINTTHISCLASS) THEN + WRITE (PUNIT,'(A,A)') 'Class ',XRDCCLASSNAMES(1) + END IF +C ------------------------------------------------------------------ +C For every residual dipolar couplig entryY... +C ------------------------------------------------------------------ + COUNT=0 + DO WHILE (COUNT.LT.RDCNUM) + COUNT=COUNT+1 +C ------------------------------------------------------------------ +C Is this the start of a new class? +C ------------------------------------------------------------------ + IF (XRDCASSNDX(CLASS).LT.COUNT) THEN + CLASS=CLASS+1 + PRINTTHISCLASS=PRINTCLASS(CLASS) + IF (XRDCASSNDX(CLASS).EQ.XRDCASSNDX(CLASS-1)) THEN + COUNT=COUNT-1 + END IF + IF (PRINTTHISCLASS) THEN + WRITE (PUNIT,'(A,A)') 'Class ',XRDCCLASSNAMES(CLASS) + END IF + END IF +C ------------------------------------------------------------------ +C If this assignment is in a class to be printed +C and make sure there is an entry for that class +C ------------------------------------------------------------------ + IF ((PRINTTHISCLASS).AND. + & (XRDCASSNDX(CLASS).NE.XRDCASSNDX(CLASS-1))) THEN +C ------------------------------------------------------------------ +C Always update RMSD +C ------------------------------------------------------------------ + CALCXRDC=XRDCCALC(COUNT) + OBSXRDC=XRDCOBS(COUNT) + ERRXRDC=XRDCERR(COUNT) + DP=(CALCXRDC-OBSXRDC) + + IF ((DP.LT.0.000).AND.(ABS(DP).GT.ERRXRDC)) THEN + DELTAXRDC=DP+ERRXRDC + ELSE IF ((DP.GT.0.000).AND.(DP.GT.ERRXRDC)) THEN + DELTAXRDC=DP-ERRXRDC + ELSE + DELTAXRDC=0.0 + END IF + + RMS=RMS+DELTAXRDC**2 + NUM=NUM+1 +C ------------------------------------------------------------------ +C Print out deviations larger than cutoff +C and update number of violations +C ------------------------------------------------------------------ + IF (ABS(DELTAXRDC).GE.CUTOFF) THEN + I=ATOMILST(ATOMIPTR(COUNT)+1) + J=ATOMJLST(ATOMJPTR(COUNT)+1) + K=ATOMKLST(ATOMKPTR(COUNT)+1) + L=ATOMLLST(ATOMLPTR(COUNT)+1) + WRITE(PUNIT,'(A,A)') '==============================', + & '==============================' + WRITE(PUNIT,'(A)') 'Set-M-atoms' + DO II=ATOMMPTR(COUNT)+1,ATOMMPTR(COUNT+1) + M=ATOMMLST(II) + WRITE(PUNIT,'(9X,4(1X,A))') SEGID(M),RESID(M), + & RES(M),TYPE(M) + END DO + WRITE(PUNIT,'(A)') 'Set-N-atoms' + DO II=ATOMNPTR(COUNT)+1,ATOMNPTR(COUNT+1) + N=ATOMNLST(II) + WRITE(PUNIT,'(9X,4(1X,A))') SEGID(N),RESID(N), + & RES(N),TYPE(N) + END DO + + WRITE(PUNIT,'(2(2X,A,1X,F8.3))') + & 'Calc: ',CALCXRDC,'Obs: ',OBSXRDC + WRITE(PUNIT,'(2X,A,1X,F8.3,2X,A,1X,F8.3)') + & 'Error: ',ERRXRDC,'Delta: ',DELTAXRDC + VIOLS=VIOLS+ONE + END IF + END IF + END DO + + IF (NUM.GT.0) THEN + RMS=SQRT(RMS / NUM) + ELSE + RMS=0.0 + ENDIF + + CALL DECLAR('RESULT','DP',' ',DUMMY2,RMS) + CALL DECLAR('VIOLATIONS','DP',' ',DUMMY2,VIOLS) + + RETURN + END diff -u -N -r xplor-nih-2.9.2/source/xdipo_rdc.fcm PARAxplor-nih-2.9.2/source/xdipo_rdc.fcm --- xplor-nih-2.9.2/source/xdipo_rdc.fcm 1970-01-01 01:00:00.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/xdipo_rdc.fcm 2003-12-05 14:34:36.000000000 +0100 @@ -0,0 +1,61 @@ +C ------------------------------------------------------------------ +C Residual dipolar coupling stuff +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + INTEGER MAXXRDCCLASSES + PARAMETER (MAXXRDCCLASSES = 40) +C ------------------------------------------------------------------ +C Arrays that hold residual dipolar coupling info +C XRDCASSNDX tells ending index of the RDC arrays for each class +C XRDCFORCES holds K1 for each class +C ------------------------------------------------------------------ + INTEGER XRDCASSNDX(MAXXRDCCLASSES) + REAL XRDCFORCES(MAXXRDCCLASSES) + REAL XRDCCOEF1(MAXXRDCCLASSES),XRDCCOEF2(MAXXRDCCLASSES) + CHARACTER*8 XRDCCLASSNAMES(MAXXRDCCLASSES) + LOGICAL PRINTCLASS(MAXXRDCCLASSES),DIPFLG +C ------------------------------------------------------------------ +C MAXXRDC is the number of slots set aside for RDC assignments +C RDCNUM is the total number of RDC entered +C ------------------------------------------------------------------ + INTEGER MAXXRDC,RDCNUM,NCLASSES,CURCLASS,XRDCMX +C ------------------------------------------------------------------ +C Pointers to arrays to hold atom numbers, observed RDC and errors +C ------------------------------------------------------------------ + INTEGER XRDCIPTR,XRDCJPTR,XRDCKPTR,XRDCLPTR,XRDCMPTR,XRDCNPTR, + & XRDCILST,XRDCJLST,XRDCKLST,XRDCLLST,XRDCMLST,XRDCNLST, + & XRDCOBSPTR,XRDCERRPTR,CALCXRDCPTR +C ------------------------------------------------------------------ +C Input modes +C ------------------------------------------------------------------ + INTEGER MODE,NEW,UPDATE + PARAMETER(NEW = 1) + PARAMETER(UPDATE = 2) +C ------------------------------------------------------------------ +C Parameters as set up in ETOR - Not used indeed +C ------------------------------------------------------------------ + DOUBLE PRECISION MCONST + PARAMETER(MCONST=0.0001D0) + DOUBLE PRECISION EPS + PARAMETER(EPS=0.1D0) +C ------------------------------------------------------------------ +C Common blocks +C ------------------------------------------------------------------ + COMMON /CXRDC/XRDCCLASSNAMES + COMMON /IXRDC1/XRDCASSNDX,MAXXRDC,RDCNUM,CURCLASS,NCLASSES,XRDCMX + COMMON /IXRDC2/XRDCIPTR,XRDCJPTR,XRDCKPTR, + & XRDCLPTR,XRDCMPTR,XRDCNPTR, + & XRDCILST,XRDCJLST,XRDCKLST, + & XRDCLLST,XRDCMLST,XRDCNLST, + & XRDCOBSPTR,XRDCERRPTR,CALCXRDCPTR,MODE + COMMON /RXRDC/XRDCFORCES,XRDCCOEF1,XRDCCOEF2 + COMMON /LXRDC/PRINTCLASS + COMMON /WXRDC/DIPFLG + + SAVE /CXRDC/ + SAVE /IXRDC1/ + SAVE /IXRDC2/ + SAVE /RXRDC/ + SAVE /LXRDC/ + SAVE /WXRDC/ diff -u -N -r xplor-nih-2.9.2/source/xfantaccr.f PARAxplor-nih-2.9.2/source/xfantaccr.f --- xplor-nih-2.9.2/source/xfantaccr.f 1970-01-01 01:00:00.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/xfantaccr.f 2003-12-05 14:34:36.000000000 +0100 @@ -0,0 +1,136 @@ +C ------------------------------------------------------------------ + SUBROUTINE FANTACCR (ATOMIPTR,ATOMJPTR,ATOMKPTR, + & ATOMILST,ATOMJLST,ATOMKLST, + & XCCROBS,XCCRERR,NCLASS,NOMEFILE) +C ------------------------------------------------------------------ +C Fits cross correlation rate coefficient on a given structure +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "psf.fcm" + INCLUDE "coord.fcm" + INCLUDE "numbers.fcm" + INCLUDE "deriv.fcm" + INCLUDE "xccr.fcm" + INCLUDE "consta.fcm" + INCLUDE "fantaxplor.fcm" + INCLUDE 'heap.fcm' + + DOUBLE PRECISION AK + LOGICAL CSWITCH,CSAVE,CERR + COMMON /CRUN/AK + COMMON /CPARAM/CSWITCH,CSAVE,CERR + INTEGER CVIOLS,IA + COMMON /CVIOL/CVIOLS + CHARACTER NOMEFILE*132 + INTEGER ATOMIPTR(*),ATOMJPTR(*),ATOMKPTR(*) + INTEGER ATOMILST(*),ATOMJLST(*),ATOMKLST(*) + + INTEGER COUNT,MM,NN,I,J,K,L,NA,M,N,CLASS,NORES + INTEGER NCLASS,NNCO,NNCC + DOUBLE PRECISION XCCROBS(*),XCCRERR(*) + DOUBLE PRECISION OBSXCCR,ERRXCCR + DOUBLE PRECISION XI,XJ,XK,YI,YJ,YK,ZI,ZJ,ZK + DOUBLE PRECISION XFA(1000),YFA(1000),ZFA(1000),XMFA(1000), + & YMFA(1000),ZMFA(1000),ERRFA(1000),OBSFA(1000), + & FENPCS(20,2000),WPROT(1000), + & XFAN(1000),YFAN(1000),ZFAN(1000) + DOUBLE PRECISION O1,O2,O3,O4,O5,O6,O7,O8,O9,O10,O11,O12 + DOUBLE PRECISION PS,DJK2,DJI2,DJK(1000),DJI(1000),ANGFACT(1000) + DOUBLE COMPLEX DUMMY2 +C ------------------------------------------------------------------ +C CLASS must reamin this during the cycle of variables exchange. +C The class is defined by NCLASS. +C ------------------------------------------------------------------ + CLASS=1 + NNCO=0 + NNCC=0 + COUNT=0 + + PRINT*,'This is FANTACROSS.' + PRINT*,'CCRNUM is:',CCRNUM + + DO WHILE(COUNT.LT.CCRNUM) + COUNT=COUNT+1 + DO WHILE ((XCCRASSNDX(CLASS).LT.COUNT).OR. + & (XCCRASSNDX(CLASS).EQ.0)) + CLASS=CLASS+1 + END DO + IF (XCCRASSNDX(CLASS).EQ.XCCRASSNDX(CLASS-1)) THEN + COUNT=COUNT-1 + ENDIF + + IF (NCLASS.EQ.CLASS) THEN + NNCO=NNCO+1 + I=ATOMIPTR(COUNT)+1 + J=ATOMJPTR(COUNT)+1 + K=ATOMKPTR(COUNT)+1 + + XI=X(ATOMILST(I)) + XJ=X(ATOMJLST(J)) + XK=X(ATOMKLST(K)) + YI=Y(ATOMILST(I)) + YJ=Y(ATOMJLST(J)) + YK=Y(ATOMKLST(K)) + ZI=Z(ATOMILST(I)) + ZJ=Z(ATOMJLST(J)) + ZK=Z(ATOMKLST(K)) + OBSFA(NNCO)=XCCROBS(COUNT) + OBSXCCR=XCCROBS(COUNT) + ERRFA(NNCO)=XCCRERR(COUNT) + ERRXCCR=XCCRERR(COUNT) + END IF + + DO MM=ATOMJPTR(COUNT)+1,ATOMJPTR(COUNT+1) + DO NN=ATOMKPTR(COUNT)+1,ATOMKPTR(COUNT+1) + IF (NCLASS.EQ.CLASS) THEN + NNCC=NNCC+1 + XMFA(NNCC)=X(ATOMILST(MM)) + XFA(NNCC)=X(ATOMJLST(NN)) + XFAN(NNCC)=X(ATOMKLST(NN)) + YMFA(NNCC)=Y(ATOMILST(MM)) + YFA(NNCC)=Y(ATOMJLST(NN)) + YFAN(NNCC)=Y(ATOMKLST(NN)) + ZMFA(NNCC)=Z(ATOMILST(MM)) + ZFA(NNCC)=Z(ATOMJLST(NN)) + ZFAN(NNCC)=Z(ATOMKLST(NN)) + NORES=ATOMJLST(NN) + NRESIDFA(NNCC)=RESID(NORES) + NRESFA(NNCC)=RES(NORES) + NTYPEFA(NNCC)=TYPE(NORES) + END IF + END DO + END DO + END DO + + IF (NNCC.EQ.NNCO) THEN + DO IA=1,CCRNUM + O1=XFA(IA)-XFAN(IA) + O2=XFA(IA)-XMFA(IA) + O3=YFA(IA)-YFAN(IA) + O4=YFA(IA)-YMFA(IA) + O5=ZFA(IA)-ZFAN(IA) + O6=ZFA(IA)-ZMFA(IA) + O7=O2**2 + O8=O4**2 + O9=O6**2 + O10=O1**2 + O11=O3**2 + O12=O5**2 + PS=O1*O2+O3*O4+O5*O6 + DJK2=O10+O11+O12 + DJK(IA)=SQRT(DJK2) + DJI2=O7+O8+O9 + DJI(IA)=SQRT(DJI2) + ANGFACT(IA)=(3.*(PS/(DJK(IA)*DJI(IA)))**2-1.) + END DO + CALL FANTALINCCR(DJI,DJK,ANGFACT,NNCO,OBSFA,ERRFA,NOMEFILE, + & WPROT,0) + ELSE + PRINT*,'Error in atom count...' + END IF + + RETURN + END diff -u -N -r xplor-nih-2.9.2/source/xfantadipo.f PARAxplor-nih-2.9.2/source/xfantadipo.f --- xplor-nih-2.9.2/source/xfantadipo.f 1970-01-01 01:00:00.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/xfantadipo.f 2003-12-05 14:34:36.000000000 +0100 @@ -0,0 +1,243 @@ +C ------------------------------------------------------------------ + SUBROUTINE FANTADIPO(ATOMIPTR,ATOMJPTR,ATOMKPTR,ATOMLPTR,ATOMMPTR, + & ATOMNPTR,ATOMILST,ATOMJLST,ATOMKLST,ATOMLLST, + & ATOMMLST,ATOMNLST,XRDCOBS,XRDCERR, + & NCLASS,NOMEFILE) +C ------------------------------------------------------------------ +C Fits residual dipolar couplings tensor on a given structure +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "psf.fcm" + INCLUDE "coord.fcm" + INCLUDE "numbers.fcm" + INCLUDE "deriv.fcm" + INCLUDE "xdipo_rdc.fcm" + INCLUDE "consta.fcm" + INCLUDE "fantaxplor.fcm" + INCLUDE 'heap.fcm' + + CHARACTER NOMEFILE*132 + DOUBLE PRECISION APARA,APERP + COMMON /FRUN/APERP,APARA + LOGICAL FDSWITCH,FDSAVE,FDERR + COMMON /FDPARAM/FDSWITCH,FDSAVE,FDERR + INTEGER FDVIOLS + COMMON /DVIOL/FDVIOLS + DOUBLE PRECISION FENERGTOT + COMMON /FENERG/FENERGTOT + DOUBLE PRECISION FDENT(20,2000),FDPER(20,2000),FDPAR(20,2000), + & FDIPOE(20,2000) + INTEGER DCONTSAVE(20) + COMMON /FDMEDIA/FDENT,FDPER,FDPAR,FDIPOE,DCONTSAVE + DOUBLE PRECISION DIPOE + COMMON /DIPOENERGY/ DIPOE + DOUBLE COMPLEX DUMMY2 + INTEGER ATOMIPTR(*),ATOMJPTR(*),ATOMKPTR(*) + INTEGER ATOMLPTR(*),ATOMMPTR(*),ATOMNPTR(*) + INTEGER ATOMILST(*),ATOMJLST(*),ATOMKLST(*) + INTEGER ATOMLLST(*),ATOMMLST(*),ATOMNLST(*) + INTEGER COUNT,MM,NN,I,J,K,L,NA,M,N,CLASS,NORES + INTEGER NCLASS,NNCO,NNCC + DOUBLE PRECISION XRDCOBS(*),XRDCERR(*) + DOUBLE PRECISION XI,XJ,XK,XL,XM,XN,YI,YJ,YK,YL,YM,YN, + & ZI,ZJ,ZK,ZL,ZM,ZN + DOUBLE PRECISION XFA(1000),YFA(1000),ZFA(1000),XMFA(1000), + & YMFA(1000),ZMFA(1000),ERRFA(1000),OBSFA(1000), + & WPROT(1000) +C ------------------------------------------------------------------ +C CLASS must remain this during the cycle of variables exchange. +C The class is defined by NCLASS. +C ------------------------------------------------------------------ + CLASS=1 + NNCO=0 + NNCC=0 + COUNT=0 + + PRINT*,'This is FANTALIN for RDC.' + PRINT*,'RDCNUM is:',RDCNUM + + DO WHILE(COUNT.LT.RDCNUM) + COUNT=COUNT+1 + DO WHILE ((XRDCASSNDX(CLASS).LT.COUNT).OR. + & (XRDCASSNDX(CLASS).EQ.0)) + CLASS=CLASS+1 + END DO + IF (XRDCASSNDX(CLASS).EQ.XRDCASSNDX(CLASS-1)) THEN + COUNT=COUNT-1 + ENDIF + + IF (NCLASS.EQ.CLASS) THEN + NNCO=NNCO+1 + I=ATOMIPTR(COUNT)+1 + J=ATOMJPTR(COUNT)+1 + K=ATOMKPTR(COUNT)+1 + L=ATOMLPTR(COUNT)+1 + XI=X(ATOMILST(I)) + XJ=X(ATOMJLST(J)) + XK=X(ATOMKLST(K)) + XL=X(ATOMLLST(L)) + YI=Y(ATOMILST(I)) + YJ=Y(ATOMJLST(J)) + YK=Y(ATOMKLST(K)) + YL=Y(ATOMLLST(L)) + ZI=Z(ATOMILST(I)) + ZJ=Z(ATOMJLST(J)) + ZK=Z(ATOMKLST(K)) + ZL=Z(ATOMLLST(L)) + OBSFA(NNCO)=XRDCOBS(COUNT) + ERRFA(NNCO)=XRDCERR(COUNT) + END IF + + DO MM=ATOMMPTR(COUNT)+1,ATOMMPTR(COUNT+1) + DO NN=ATOMNPTR(COUNT)+1,ATOMNPTR(COUNT+1) + IF (NCLASS.EQ.CLASS) THEN + NNCC=NNCC+1 + XMFA(NNCC)=X(ATOMMLST(MM)) + XFA(NNCC)=X(ATOMNLST(NN)) + YMFA(NNCC)=Y(ATOMMLST(MM)) + YFA(NNCC)=Y(ATOMNLST(NN)) + ZMFA(NNCC)=Z(ATOMMLST(MM)) + ZFA(NNCC)=Z(ATOMNLST(NN)) + NORES=ATOMNLST(NN) + NRESIDFA(NNCC)=RESID(NORES) + NRESFA(NNCC)=RES(NORES) + NTYPEFA(NNCC)=TYPE(NORES) + END IF + END DO + END DO + END DO + + IF (NNCC.EQ.NNCO) THEN + CALL FANTALINEXE(XFA,YFA,ZFA,XMFA,YMFA,ZMFA,NNCO,OBSFA,ERRFA, + & 0,NOMEFILE,WPROT,0,DIPFLG) + PRINT*,'FDERR is=',FDERR + IF (FDERR.EQ.1) THEN + CALL FANTAERRORE(XFA,YFA,ZFA,XMFA,YMFA,ZMFA,NNCO,OBSFA, + & ERRFA,0,DIPFLG) + END IF + ELSE + PRINT*,'Error in atom count...' + END IF + + PRINT*,'FDSWITCH is:',FDSWITCH + IF (FDSWITCH.EQ.1) THEN + XRDCCOEF1(NCLASS)=APARA + XRDCCOEF2(NCLASS)=APERP + PRINT*,'For class:',NCLASS + PRINT*,'New value for A1 =',APARA + PRINT*,'New value for A2 =',APERP + END IF + + IF (FDSAVE.EQ.1) THEN + DCONTSAVE(NCLASS)=DCONTSAVE(NCLASS)+1 + FDENT(NCLASS,DCONTSAVE(NCLASS))=FENERGTOT + FDIPOE(NCLASS,DCONTSAVE(NCLASS))=DIPOE + FDPAR(NCLASS,DCONTSAVE(NCLASS))=APARA + FDPER(NCLASS,DCONTSAVE(NCLASS))=APERP + PRINT*,'XDIPO_RDC energy =',DIPOE + PRINT*,'Total energy =',FENERGTOT + PRINT*,'DCONTSAVE is:',DCONTSAVE(NCLASS) + END IF + + CALL DECLAR('DCHIAX','DP',' ', DUMMY2,APARA) + CALL DECLAR('DCHIRH','DP',' ', DUMMY2,APERP) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE FANTADIPOMEDIA (FMA,PERC) +C ------------------------------------------------------------------ +C Averages tensor parameters on a given sub-ensemble of structures +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + DOUBLE PRECISION VETE(2,2000),VETT(3,2000) + DOUBLE PRECISION FDENT(20,2000),FDPER(20,2000),FDPAR(20,2000), + & FDIPOE(20,2000) + COMMON /FDMEDIA/ FDENT,FDPER,FDPAR,FDIPOE,DCONTSAVE + INTEGER DCONTSAVE(20) + INTEGER FMA,COUNT,CB,NUM,COUNT1,CONFPRE + DOUBLE PRECISION PERC + DOUBLE PRECISION M1PAR,M1PER,DEVPAR + DOUBLE COMPLEX DUMMY2 +C ------------------------------------------------------------------ +C Through FMA the class on which the average must be done is passed. +C If it is zero, all the counters are reset. +C ------------------------------------------------------------------ + IF (FMA.EQ.0) THEN + DO COUNT=1,19 + DCONTSAVE(COUNT)=0 + END DO + PRINT*,'Mean vectors reset.' + ELSE +C ------------------------------------------------------------------ +C The best are selected +C ------------------------------------------------------------------ + DO COUNT=1,DCONTSAVE(FMA) + VETE(1,COUNT)=0 + END DO + + DO COUNT=1,DCONTSAVE(FMA) + NUM=DCONTSAVE(FMA) + DO WHILE (VETE(1,NUM).EQ.1) + NUM=NUM-1 + END DO + DO COUNT1=1,DCONTSAVE(FMA)-1 + IF ((FDENT(FMA,COUNT1).LE.FDENT(FMA,NUM)).AND. + & (VETE(1,COUNT1).EQ.0.0)) THEN + NUM=COUNT1 + END IF + END DO + PRINT*,'NUM is:',NUM + VETE(1,NUM)=1 + VETE(2,COUNT)=FDENT(FMA,NUM) + VETT(1,COUNT)=FDPAR(FMA,NUM) + VETT(2,COUNT)=FDPER(FMA,NUM) + END DO + + CONFPRE=NINT((DBLE(DCONTSAVE(FMA))*(PERC/100.))) + IF (CONFPRE.EQ.0) THEN + CONFPRE=1 + END IF + PRINT*,'CONFPRE is:',CONFPRE + PRINT*,'PERC is:',PERC + PRINT*,'DBLE(CONTSAVE(FMA)) is:',DBLE(DCONTSAVE(FMA)) + + DO COUNT=1,DCONTSAVE(FMA) + PRINT*,'Number',COUNT,'Energy',VETE(1,COUNT),VETE(2,COUNT) + PRINT*,'A1',VETT(1,COUNT),' A2',VETT(2,COUNT) + END DO + + CB=0 + M1PAR=0 + M1PER=0 + DO COUNT=1,CONFPRE + M1PAR=M1PAR+VETT(1,COUNT) + M1PER=M1PER+VETT(2,COUNT) + CB=CB+1 + PRINT*,'A1 =',VETT(1,COUNT),'A2 =',VETT(2,COUNT) + END DO + M1PAR=M1PAR/CB + M1PER=M1PER/CB + + CB=0 + DEVPAR=0 + DO COUNT=1,CONFPRE + DEVPAR=DEVPAR+(VETT(1,COUNT)-M1PAR)**2 + END DO + DEVPAR=DEVPAR/CONFPRE + PRINT*,'Averaged new A1 =',M1PAR + PRINT*,'Averaged new A2 =',M1PER + PRINT*,'Standard deviation for new A1 =',SQRT(DEVPAR) + + CALL DECLAR('DCHIAX','DP',' ',DUMMY2,M1PAR) + CALL DECLAR('DCHIRH','DP',' ',DUMMY2,M1PER) + + END IF + + RETURN + END diff -u -N -r xplor-nih-2.9.2/source/xplorFunc.f PARAxplor-nih-2.9.2/source/xplorFunc.f --- xplor-nih-2.9.2/source/xplorFunc.f 2003-07-22 14:58:51.000000000 +0200 +++ PARAxplor-nih-2.9.2/source/xplorFunc.f 2003-12-05 14:34:36.000000000 +0100 @@ -91,6 +91,26 @@ C initialize Dipolar coupling refinement C CALL XDIPINIT +c ------------------------------------------------------------------ +c PARArestraints (2003) +c +C +C initialize Pseudocontact shifts refinement +C + CALL XPCSINIT +C +C initialize Residual dipolar couplings refinement +C + CALL XRDCINIT +C +C initialize Cross-correlation rates refinement +C + CALL XCCRINIT +C +C initialize RDC-ANGLES refinement +C + CALL XANGLESINIT +c ------------------------------------------------------------------ C C initialize Diffusion Anisotropy refinement C @@ -327,6 +347,15 @@ &' SANIsotropy END ', &' DANIsotropy END ', &' DIPOlar END ', +c ------------------------------------------------------------------ +c Commands related to PARArestraints (2003) +c + &' XPCShift END ', + &' XRDCouplings END ', + &' XCCRates END ', + &' XANGles END ', + &' XT1Distances END ', +c ------------------------------------------------------------------ &' DCSAnisotropy END ', &' HBDA END ', &' VEAN END ', @@ -368,6 +397,20 @@ CALL READSANI ELSE IF (WD(1:4).EQ.'DIPO' .or. WD(1:4).EQ.'XDIP' ) THEN CALL READXDIP +c ------------------------------------------------------------------ +c Commands for PARArestraints (2003) +c + ELSE IF (WD(1:4).EQ.'XPCS' ) THEN + CALL READXPCS + ELSE IF (WD(1:4).EQ.'XRDC' ) THEN + CALL READXRDC + ELSE IF (WD(1:4).EQ.'XCCR' ) THEN + CALL READXCCR + ELSE IF (WD(1:4).EQ.'XANG' ) THEN + CALL READXANGLES + ELSE IF (WD(1:4).EQ.'XT1D' ) THEN + CALL XT1DIST +c ------------------------------------------------------------------ ELSE IF (WD(1:4).EQ.'DCSA' ) THEN CALL READDCSA ELSE IF (WD(1:4).EQ.'HBDA' ) THEN diff -u -N -r xplor-nih-2.9.2/source/xt1dist.f PARAxplor-nih-2.9.2/source/xt1dist.f --- xplor-nih-2.9.2/source/xt1dist.f 1970-01-01 01:00:00.000000000 +0100 +++ PARAxplor-nih-2.9.2/source/xt1dist.f 2003-12-05 14:34:36.000000000 +0100 @@ -0,0 +1,215 @@ +C ------------------------------------------------------------------ + SUBROUTINE XT1DIST +C ------------------------------------------------------------------ +C Reads in relaxation rate enhancement information +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + IMPLICIT NONE + + INCLUDE "comand.fcm" + INCLUDE "funct.fcm" + INCLUDE "psf.fcm" + INCLUDE "heap.fcm" + INCLUDE "numbers.fcm" + + INTEGER SPTR + CHARACTER KT1INPUT*132,KT1OUTPUT*132,KPDB*132 + DOUBLE PRECISION KT1 + + SPTR=ALLHP(INTEG4(NATOM)) + CALL PUSEND('XT1D>') +862 CONTINUE + CALL NEXTWD('XT1D>') + CALL MISCOM('XT1D>',USED) + IF (.NOT.USED) THEN +C ------------------------------------------------------------------ +C Documentation +C ------------------------------------------------------------------ + IF (WD(1:4).EQ.'HELP') THEN + WRITE(DUNIT,'(10X,A)') + & ' XT1D {} END ', + & ' :== ', + & ' KINPut ', + & ' * InputFile Kvalue OutputFile *', + & ' STRUcture ', + & ' * PDBfile InputFile OutputFile *' +C ------------------------------------------------------------------ +C Restraint conversion using a given value of K +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'KINP') THEN + CALL NEXTFI('Input file for T1 =',KT1INPUT) + PRINT*,'Using as input file for T1:',KT1INPUT + CALL NEXTF('Value of K =',KT1) + PRINT*,'Using as value of K:',KT1 + CALL NEXTFI('Output NOE-like file for T1 =',KT1OUTPUT) + PRINT*,'Using as output NOE-like file for T1:',KT1OUTPUT + CALL KT1TYPECALC(KT1INPUT,KT1,KT1OUTPUT) +C ------------------------------------------------------------------ +C Restraint conversion using a structure +C ------------------------------------------------------------------ + ELSE IF (WD(1:4).EQ.'STRU') THEN + CALL NEXTFI('Input PDB file for T1 =',KPDB) + PRINT*,'Using as input PDB file for T1:',KPDB + CALL NEXTFI('Input file for T1 =',KT1INPUT) + PRINT*,'Using as input file for T1:',KT1INPUT + CALL NEXTFI('Output NOE-like file for T1 =',KT1OUTPUT) + PRINT*,'Using as output NOE-like file for T1:',KT1OUTPUT + CALL KT1PDB(KT1INPUT,KPDB,KT1OUTPUT) +C ------------------------------------------------------------------ +C Check for END statement +C ------------------------------------------------------------------ + ELSE + CALL CHKEND('XT1D>',DONE) + END IF + END IF + IF (.NOT.(DONE)) GOTO 862 + DONE=.FALSE. + CALL FREHP(SPTR,INTEG4(NATOM)) + + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE KT1TYPECALC (KT1INPUT,KT1,KT1OUTPUT) +C ------------------------------------------------------------------ +C Converts T1 data to distance restraints using a given value of K +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + CHARACTER KT1INPUT*132,KT1OUTPUT*132 + CHARACTER BEGIN*3 + DOUBLE PRECISION KT1,R,RNOE + INTEGER N1A,N2A + CHARACTER*4 R1A,R2A + INTEGER NLINES + + OPEN(21,FILE=KT1INPUT,STATUS='OLD',ERR=666) + OPEN(31,FILE=KT1OUTPUT) + + NLINES=0 + DO I=1,10000 +10 READ(21,91,END=20,ERR=21) BEGIN + NLINES=NLINES+1 + IF (BEGIN.NE.' ') THEN + IF (BEGIN.EQ.'KT1') THEN + BACKSPACE 21 + READ(21,92) BEGIN,N1A,R1A,N2A,R2A,R + RNOE=((KT1/R)**(1/6.))+1 + WRITE(31,94) ' assign (resid ',N1A,' and name ',R1A, + & ') (resid ',N2A,' and name ',R2A,') ', + & RNOE,RNOE,' 0' + END IF + END IF +20 CONTINUE + END DO +21 CONTINUE + CLOSE (21) + CLOSE (31) + +91 FORMAT (A3) +92 FORMAT (A3,1X,I4,1X,A4,1X,I4,1X,A4,1X,F6.2) +94 FORMAT (A15,I4,A10,A4,A9,I4,A10,A4,A2,F5.1,F5.1,A2) + + RETURN +666 PRINT*,'The file:',KT1INPUT,'does not exist...' + RETURN + END +C ------------------------------------------------------------------ + SUBROUTINE KT1PDB (KT1INPUT,KPDB,KT1OUTPUT) +C ------------------------------------------------------------------ +C Converts T1 data to distance restraints using a given structure +C +C By Gabriele Cavallaro, Andrea Giachetti and Giacomo Parigi (2003) +C ------------------------------------------------------------------ + CHARACTER KT1INPUT*132,KT1OUTPUT*132,KPDB*132 + CHARACTER BEGIN*4,BEGIN1*6 + DOUBLE PRECISION KT1,R(2000),RNOE,RMAX + DOUBLE PRECISION X1,Y1,Z1,RA,X2,Y2,Z2 + INTEGER N1A(2000),N2A(2000) + CHARACTER*4 R1A(2000),R2A(2000) + INTEGER NLINES,NUMBERATOM,NORES,BB + CHARACTER NAMEATOM*4 + + OPEN(21,FILE=KT1INPUT,STATUS='OLD',ERR=666) + + RMAX=0 + NLINES=0 + DO I=1,10000 +10 READ(21,91,END=20,ERR=21) BEGIN + IF (BEGIN.NE.' ') THEN + IF (BEGIN.EQ.'KT1 ') THEN + NLINES=NLINES+1 + BACKSPACE 21 + READ(21,92) BEGIN,N1A(NLINES),R1A(NLINES),N2A(NLINES), + & R2A(NLINES),R(NLINES) + END IF + END IF +20 CONTINUE + END DO +21 CONTINUE + CLOSE(21) + + K=0 + DO J=1,NLINES + OPEN(11,FILE=KPDB,STATUS='OLD',ERR=667) + BB=0 + DO I=1,10000 + READ(11,91,END=102) BEGIN1 + IF (BEGIN1.NE.' ') THEN + IF (BEGIN1.EQ.'ATOM ') THEN + BACKSPACE 11 + READ(11,'(5X,I6,1X,A4,6X,I4,4X,3F8.3)') + & NUMBERATOM,NAMEATOM,NORES,XA,YA,ZA + IF (NORES.EQ.N1A(J)) THEN + IF (NAMEATOM.EQ.R1A(J)) THEN + X1=XA + Y1=YA + Z1=ZA + BB=BB+1 + END IF + END IF + IF (NORES.EQ.N2A(J)) THEN + IF (NAMEATOM.EQ.R2A(J)) THEN + X2=XA + Y2=YA + Z2=ZA + BB=BB+1 + END IF + END IF + END IF + END IF + END DO +102 CONTINUE + IF (BB.NE.2) THEN + PRINT*,'%T1DIST-ERR: Mismatch in PDB file at line:',J + STOP + END IF + RA=SQRT((X1-X2)**2+(Y1-Y2)**2+(Z1-Z2)**2) + K=(RA**6)*R(J) + IF (K.GT.KMAX) THEN + KMAX=K + END IF + CLOSE (11) + END DO + PRINT*,'Value of KMAX:',KMAX + + OPEN(31,FILE=KT1OUTPUT) + DO J=1,NLINES + RNOE=(KMAX/R(J))**(1/6.) + WRITE(31,94) ' assign (resid ',N1A(J),' and name ',R1A(J), + & ') (resid ',N2A(J),' and name ',R2A(J),') ', + & RNOE,RNOE,' 0' + END DO + CLOSE(31) + +91 FORMAT (A4) +92 FORMAT (A3,1X,I4,1X,A4,1X,I4,1X,A4,1X,F6.2) +94 FORMAT (A15,I4,A10,A4,A9,I4,A10,A4,A2,F8.3,F8.3,A2) + + RETURN +666 PRINT*,'The file:',KT1INPUT,' does not exist...' + RETURN +667 PRINT*,'The file:',KPDB,' does not exist...' + RETURN + + END