From: Cezary Czaplewski Date: Thu, 24 Mar 2016 02:08:05 +0000 (+0100) Subject: cluster_wham src-M agrees with src for single chain X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=3b4fa21bffe34c0d670a637e96a9b275f8d04ac2 cluster_wham src-M agrees with src for single chain --- diff --git a/source/cluster/wham/src-M/COMMON.CLUSTER b/source/cluster/wham/src-M/COMMON.CLUSTER index 4477d19..f1ad0fd 100644 --- a/source/cluster/wham/src-M/COMMON.CLUSTER +++ b/source/cluster/wham/src-M/COMMON.CLUSTER @@ -4,11 +4,11 @@ real*4 diss,allcart double precision enetb,entfac,totfree,energy,rmstb integer ncut,ngr,licz,nconf,iass,icc,mult,list_conf, - & nss_all,ihpb_all,jhpb_all,iass_tot,iscore,nprop + & nss_all,ihpb_all,jhpb_all,iass_tot,iscore,nprop,nclust common /clu/ diss(maxdist),energy(0:maxconf), & enetb(max_ene,maxstr_proc),ecut, & entfac(maxconf),totfree(0:maxconf),totfree_gr(maxgr), - & rcutoff(max_cut+1),ncut,min_var,tree,plot_tree,lgrp + & rcutoff(max_cut+1),ncut,nclust,min_var,tree,plot_tree,lgrp common /clu1/ ngr,licz(maxgr),nconf(maxgr,maxingr),iass(maxgr), & iass_tot(maxgr,max_cut),list_conf(maxconf) common /alles/ allcart(3,maxres2,maxstr_proc),rmstb(maxconf), diff --git a/source/cluster/wham/src-M/energy_p_new.F b/source/cluster/wham/src-M/energy_p_new.F index c7c2c2f..4551e68 100644 --- a/source/cluster/wham/src-M/energy_p_new.F +++ b/source/cluster/wham/src-M/energy_p_new.F @@ -3349,7 +3349,8 @@ C etheta=0.0D0 c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) do i=ithet_start,ithet_end - if (itype(i-1).eq.21) cycle + if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. + &(itype(i).eq.ntyp1)) cycle dethetai=0.0d0 dephii=0.0d0 dephii1=0.0d0 @@ -3359,7 +3360,7 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3 .and. itype(i-2).ne.21) then + if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -3373,13 +3374,13 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) enddo else phii=0.0d0 - ityp1=nthetyp+1 + ityp1=ithetyp(itype(i-2)) do k=1,nsingle cosph1(k)=0.0d0 sinph1(k)=0.0d0 enddo endif - if (i.lt.nres .and. itype(i).ne.21) then + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -3394,7 +3395,7 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) enddo else phii1=0.0d0 - ityp3=nthetyp+1 + ityp3=ithetyp(itype(i)) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0 diff --git a/source/cluster/wham/src-M/main_clust.F b/source/cluster/wham/src-M/main_clust.F index a2e4769..d33ce25 100644 --- a/source/cluster/wham/src-M/main_clust.F +++ b/source/cluster/wham/src-M/main_clust.F @@ -23,7 +23,7 @@ C logical printang(max_cut) integer printpdb(max_cut) integer printmol2(max_cut) - character*240 lineh + character*240 lineh,scrachdir2d REAL CRIT(maxconf),MEMBR(maxconf) REAL CRITVAL(maxconf-1) INTEGER IA(maxconf),IB(maxconf) @@ -34,13 +34,14 @@ C DIMENSION NN(maxconf),DISNN(maxconf) LOGICAL FLAG(maxconf) integer i,j,k,l,m,n,len,lev,idum,ii,ind,ioffset,jj,icut,ncon, - & it,ncon_work,ind1 + & it,ncon_work,ind1,ilen,is,ie double precision t1,t2,tcpu,difconf real diss_(maxdist) double precision varia(maxvar) double precision hrtime,mintime,sectime logical eof + external ilen #ifdef MPI call MPI_Init( IERROR ) call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR ) @@ -80,6 +81,12 @@ c write (iout,*) "after permut" c call flush(iout) print *,'MAIN: nnt=',nnt,' nct=',nct + if (nclust.gt.0) then + PRINTANG(1)=.TRUE. + PRINTPDB(1)=outpdb + printmol2(1)=outmol2 + ncut=0 + else DO I=1,NCUT PRINTANG(I)=.FALSE. PRINTPDB(I)=0 @@ -91,12 +98,21 @@ c call flush(iout) printmol2(i)=outmol2 ENDIF ENDDO + endif + if (ncut.gt.0) then write (iout,*) 'Number of cutoffs:',NCUT write (iout,*) 'Cutoff values:' DO ICUT=1,NCUT WRITE(IOUT,'(8HRCUTOFF(,I2,2H)=,F8.1,2i2)')ICUT,RCUTOFF(ICUT), & printpdb(icut),printmol2(icut) ENDDO + else if (nclust.gt.0) then + write (iout,'("Number of clusters requested",i5)') nclust + else + if (me.eq.Master) + & write (iout,*) "ERROR: Either nclust or ncut must be >0" + stop + endif DO I=1,NRES-3 MULT(I)=1 ENDDO @@ -112,7 +128,6 @@ c call flush(iout) #ifdef MPI call work_partition(.true.,ncon) #endif - call probabl(iT,ncon_work,ncon,*20) if (ncon_work.lt.2) then @@ -174,7 +189,8 @@ c write (iout,'(2i4,i10,f10.5)') i,j,ind,DISS(IND) & scount(0),idispl(0),MPI_REAL,Master,MPI_COMM_WORLD, IERROR) if (me.eq.master) then #endif - open(80,file='/tmp/distance',form='unformatted') + scrachdir2d=scratchdir(:ilen(scratchdir))//'distance' + open(80,file=scrachdir2d,form='unformatted') do i=1,ndis write(80) diss(i) enddo @@ -247,29 +263,39 @@ C CALL HCASS(N,IA,IB,CRIT,LEV,ICLASS,HVALS,IORDER,CRITVAL,HEIGHT) c CALL HCDEN(LEV,IORDER,HEIGHT,CRITVAL) +c 3/3/16 AL: added explicit number of cluters + if (nclust.gt.0) then + is=nclust-1 + ie=nclust-1 + icut=1 + else + is=1 + ie=lev-1 + endif do i=1,maxgr licz(i)=0 enddo icut=1 - i=1 - NGR=i+1 + i=is + NGR=is+1 do j=1,n licz(iclass(j,i))=licz(iclass(j,i))+1 nconf(iclass(j,i),licz(iclass(j,i)))=j c write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)), c & nconf(iclass(j,i),licz(iclass(j,i))) enddo - do i=1,lev-1 - +c do i=1,lev-1 + do i=is,ie idum=lev-i DO L=1,LEV IF (HEIGHT(L).EQ.IDUM) GOTO 190 ENDDO 190 IDUM=L - write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM), - & " icut",icut," cutoff",rcutoff(icut) - IF (CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN - WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut) +c write(IOUT,*) "i+1",i+1," idum",idum," critval",CRITVAL(IDUM), +c & " icut",icut," cutoff",rcutoff(icut) + IF (nclust.gt.0.or.CRITVAL(IDUM).LT. RCUTOFF(ICUT)) THEN + if (nclust.le.0) + & WRITE (iout,'(/a,f10.5)') 'AT CUTOFF:',rcutoff(icut) write (iout,'(a,f8.2)') 'Maximum distance found:', & CRITVAL(IDUM) CALL SRTCLUST(ICUT,ncon_work,iT) @@ -282,9 +308,10 @@ c & nconf(iclass(j,i),licz(iclass(j,i))) do l=1,maxgr licz(l)=0 enddo + ii=i-is+1 do j=1,n - licz(iclass(j,i))=licz(iclass(j,i))+1 - nconf(iclass(j,i),licz(iclass(j,i)))=j + licz(iclass(j,ii))=licz(iclass(j,ii))+1 + nconf(iclass(j,ii),licz(iclass(j,ii)))=j c write (iout,*) i,j,iclass(j,i),licz(iclass(j,i)), c & nconf(iclass(j,i),licz(iclass(j,i))) cd print *,j,iclass(j,i), @@ -310,17 +337,17 @@ C C close(icbase,status="delete") #ifdef MPI - call MPI_Finalize(MPI_COMM_WORLD,IERROR) + call MPI_Finalize(IERROR) #endif stop '********** Program terminated normally.' 20 write (iout,*) "Error reading coordinates" #ifdef MPI - call MPI_Finalize(MPI_COMM_WORLD,IERROR) + call MPI_Finalize(IERROR) #endif stop 30 write (iout,*) "Error reading reference structure" #ifdef MPI - call MPI_Finalize(MPI_COMM_WORLD,IERROR) + call MPI_Finalize(IERROR) #endif stop end diff --git a/source/cluster/wham/src-M/read_coords.F b/source/cluster/wham/src-M/read_coords.F index c34aca4..f1e2c7b 100644 --- a/source/cluster/wham/src-M/read_coords.F +++ b/source/cluster/wham/src-M/read_coords.F @@ -212,6 +212,7 @@ c call flush(iout) enddo enddo else + itmp=0 #if (defined(AIX) && !defined(JUBL)) call xdrf3dfcoord_(ixdrf, csingle, itmp, prec, iret) if (iret.eq.0) goto 101 diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F index 4c15495..e33b6d6 100644 --- a/source/cluster/wham/src-M/readrtns.F +++ b/source/cluster/wham/src-M/readrtns.F @@ -38,7 +38,8 @@ C min_var=(index(controlcard,'MINVAR').gt.0) plot_tree=(index(controlcard,'PLOT_TREE').gt.0) punch_dist=(index(controlcard,'PUNCH_DIST').gt.0) - call readi(controlcard,'NCUT',ncut,1) + call readi(controlcard,'NCUT',ncut,0) + call readi(controlcard,'NCLUST',nclust,5) call readi(controlcard,'SYM',symetr,1) write (iout,*) 'sym', symetr call readi(controlcard,'NSTART',nstart,0) @@ -49,7 +50,8 @@ C lgrp=(index(controlcard,'LGRP').gt.0) caonly=(index(controlcard,'CA_ONLY').gt.0) print_dist=(index(controlcard,'PRINT_DIST').gt.0) - call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0) + if (ncut.gt.0) + & call multreada(controlcard,'CUTOFF',rcutoff,ncut,-1.0d0) call readi(controlcard,'IOPT',iopt,2) lside = index(controlcard,"SIDE").gt.0 efree = index(controlcard,"EFREE").gt.0 @@ -206,6 +208,19 @@ C Convert sequence to numeric code do i=1,nres itype(i)=rescode(i,sequence(i),iscode) enddo + if (itype(2).eq.10) then + write (iout,*) + & "Glycine is the first full residue, initial dummy deleted" + do i=1,nres + itype(i)=itype(i+1) + enddo + nres=nres-1 + endif + if (itype(nres).eq.10) then + write (iout,*) + & "Glycine is the last full residue, terminal dummy deleted" + nres=nres-1 + endif print *,nres print '(20i4)',(itype(i),i=1,nres) diff --git a/source/cluster/wham/src-M/wrtclust.f b/source/cluster/wham/src-M/wrtclust.f index f2f3eb7..3915ebc 100644 --- a/source/cluster/wham/src-M/wrtclust.f +++ b/source/cluster/wham/src-M/wrtclust.f @@ -84,12 +84,16 @@ C 12/8/93 Estimation of "diameters" of the subsequent families. ave_dim=0.0 amax_dim=0.0 c write (iout,*) "ecut",ecut + emin=totfree(nconf(igr,1)) +c write (2,*) "emin",emin," ecut",ecut do i=2,licz(igr) ii=nconf(igr,i) +c write (2,*) " igr",igr," i",i," ii",ii," totfree",totfree(ii), +c & " emin",emin," diff",totfree(ii)-emin," ecut",ecut if (totfree(ii)-emin .gt. ecut) goto 10 do j=1,i-1 jj=nconf(igr,j) - if (jj.eq.1) exit +c if (jj.eq.1) exit if (ii.lt.jj) then ind=ioffset(ncon,ii,jj) else @@ -112,9 +116,12 @@ c & list_conf(jj),curr_dist & '; average distance in the family:',ave_dim rmsave(igr)=0.0d0 qpart=0.0d0 + emin=totfree(nconf(igr,1)) do i=1,licz(igr) icon=nconf(igr,i) - boltz=dexp(-totfree(icon)) + boltz=dexp(-totfree(icon)+emin) +c write (2,*) "igr",igr," i",i," icon",icon," totfree", +c & totfree(icon)," emin",emin," boltz",boltz," rms",rmstb(icon) rmsave(igr)=rmsave(igr)+boltz*rmstb(icon) qpart=qpart+boltz enddo @@ -191,7 +198,7 @@ c Write conformations of the family i to PDB files c write (iout,*) i,ncon_out,nconf(i,ncon_out), c & totfree(nconf(i,ncon_out)),emin1,ecut enddo - write (iout,*) "ncon_out",ncon_out +c write (iout,*) "ncon_out",ncon_out call flush(iout) do j=1,nres tempfac(1,j)=5.0d0 diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index af921d0..709a471 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -3399,7 +3399,8 @@ C etheta=0.0D0 c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) do i=ithet_start,ithet_end - if (itype(i-1).eq.21) cycle + if ((itype(i-1).eq.ntyp1).or.(itype(i-2).eq.ntyp1).or. + &(itype(i).eq.ntyp1)) cycle dethetai=0.0d0 dephii=0.0d0 dephii1=0.0d0 @@ -3409,7 +3410,7 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) coskt(k)=dcos(k*theti2) sinkt(k)=dsin(k*theti2) enddo - if (i.gt.3 .and. itype(i-2).ne.21) then + if (i.gt.3 .and. itype(max0(i-3,1)).ne.ntyp1) then #ifdef OSF phii=phi(i) if (phii.ne.phii) phii=150.0 @@ -3423,13 +3424,13 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) enddo else phii=0.0d0 - ityp1=nthetyp+1 + ityp1=ithetyp(itype(i-2)) do k=1,nsingle cosph1(k)=0.0d0 sinph1(k)=0.0d0 enddo endif - if (i.lt.nres .and. itype(i).ne.21) then + if (i.lt.nres .and. itype(i+1).ne.ntyp1) then #ifdef OSF phii1=phi(i+1) if (phii1.ne.phii1) phii1=150.0 @@ -3444,7 +3445,7 @@ c write (iout,*) "ithetyp",(ithetyp(i),i=1,ntyp1) enddo else phii1=0.0d0 - ityp3=nthetyp+1 + ityp3=ithetyp(itype(i)) do k=1,nsingle cosph2(k)=0.0d0 sinph2(k)=0.0d0