From 549a47d7bb5ff7b75c531eec0c406ae46281e1bb Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Fri, 4 Mar 2016 17:57:37 +0900 Subject: [PATCH] Adam's changes to wham and cluster following previous commit XG.. and ..GX sequence is automatically corrected also in cluster_wham it is possible to select number of clusters instead of cutoff Conflicts: source/cluster/wham/src/COMMON.CONTROL source/cluster/wham/src/energy_p_new.F source/cluster/wham/src/probabl.F source/cluster/wham/src/readpdb.F source/cluster/wham/src/readrtns.F source/wham/src/energy_p_new.F --- source/cluster/wham/src/COMMON.CLUSTER | 4 +- source/cluster/wham/src/COMMON.CONTROL | 14 + source/cluster/wham/src/energy_p_new.F | 6 + source/cluster/wham/src/main_clust.F | 48 +++- source/cluster/wham/src/probabl.F | 14 +- source/cluster/wham/src/readrtns.F | 489 +++++++++++++++++++++++++++++++- source/cluster/wham/src/wrtclust.f | 11 +- source/wham/src/energy_p_new.F | 6 + source/wham/src/molread_zs.F | 13 + 9 files changed, 583 insertions(+), 22 deletions(-) diff --git a/source/cluster/wham/src/COMMON.CLUSTER b/source/cluster/wham/src/COMMON.CLUSTER index 4477d19..f1ad0fd 100644 --- a/source/cluster/wham/src/COMMON.CLUSTER +++ b/source/cluster/wham/src/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/COMMON.CONTROL b/source/cluster/wham/src/COMMON.CONTROL index 8c9e317..c195fff 100644 --- a/source/cluster/wham/src/COMMON.CONTROL +++ b/source/cluster/wham/src/COMMON.CONTROL @@ -2,8 +2,22 @@ integer iscode,indpdb,outpdb,outmol2,iopt,nstart,nend,constr_dist logical refstr,pdbref,punch_dist,print_dist,caonly,lside, & lprint_cart,lprint_int,from_cart,efree,from_bx,from_cx, +<<<<<<< HEAD & with_dihed_constr common /cntrl/ betaT,iscode,indpdb,refstr,pdbref,outpdb,outmol2, & punch_dist,print_dist,caonly,lside,lprint_cart,lprint_int, & from_cart,from_bx,from_cx,efree,iopt,nstart,nend,constr_dist, & with_dihed_constr +======= + & with_dihed_constr,out1file,print_homology_restraints, + & print_contact_map,print_homology_models + common /cntrl/ betaT,iscode,indpdb,refstr,pdbref,outpdb,outmol2, + & punch_dist,print_dist,caonly,lside,lprint_cart,lprint_int, + & from_cart,from_bx,from_cx,efree,iopt,nstart,nend,constr_dist, + & with_dihed_constr, constr_homology,homol_nset,out1file, + & print_contact_map + common /cntrlr/ waga_homology(10), + & waga_dist, waga_angle, waga_theta, waga_d, dist_cut,iset, + & print_homology_restraints,print_homology_models + +>>>>>>> 3d6f9e7... Adam's changes to wham and cluster following previous commit diff --git a/source/cluster/wham/src/energy_p_new.F b/source/cluster/wham/src/energy_p_new.F index 636f983..20c43f5 100644 --- a/source/cluster/wham/src/energy_p_new.F +++ b/source/cluster/wham/src/energy_p_new.F @@ -4509,6 +4509,12 @@ c 3 = SC...Ca...Ca...SCi cosphi=dcos(j*tauangle(intertyp,i)) sinphi=dsin(j*tauangle(intertyp,i)) esccor=esccor+v1ij*cosphi+v2ij*sinphi +<<<<<<< HEAD +======= +#ifdef DEBUG + esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi +#endif +>>>>>>> 3d6f9e7... Adam's changes to wham and cluster following previous commit gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) enddo gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci diff --git a/source/cluster/wham/src/main_clust.F b/source/cluster/wham/src/main_clust.F index 4fe3789..5a61e8d 100644 --- a/source/cluster/wham/src/main_clust.F +++ b/source/cluster/wham/src/main_clust.F @@ -34,7 +34,7 @@ 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,ilen + & it,ncon_work,ind1,ilen,is,ie double precision t1,t2,tcpu,difconf real diss_(maxdist) @@ -76,6 +76,12 @@ c if (refstr) call read_ref_structure(*30) 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 @@ -87,12 +93,21 @@ c if (refstr) call read_ref_structure(*30) 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 @@ -243,29 +258,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) @@ -278,9 +303,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), diff --git a/source/cluster/wham/src/probabl.F b/source/cluster/wham/src/probabl.F index 53735e8..fa10f58 100644 --- a/source/cluster/wham/src/probabl.F +++ b/source/cluster/wham/src/probabl.F @@ -176,16 +176,20 @@ cc if (wcorr6.eq.0) ecorr6=0.0d0 #ifdef MPI Fdimless_(i)=beta_h(ib)*etot+entfac(ii) totfree_(i)=etot +#ifdef DEBUG + write (iout,*) "etrop", i,ii,ib, + & 1.0d0/(1.987d-3*beta_h(ib)),totfree(i), + & entfac(ii),Fdimless_(i) +#endif #else Fdimless(i)=beta_h(ib)*etot+entfac(ii) totfree(i)=etot -#endif #ifdef DEBUG - write (iout,*) "etrop", i,ii,ib, & 1.0d0/(1.987d-3*beta_h(ib)),totfree(i), & entfac(ii),Fdimless(i) #endif +#endif enddo ! i #ifdef MPI call MPI_Gatherv(Fdimless_(1),scount(me), @@ -205,7 +209,11 @@ cc if (wcorr6.eq.0) ecorr6=0.0d0 entfac(i)=entfac_(i) enddo #endif +<<<<<<< HEAD c#define DEBUG +======= + +>>>>>>> 3d6f9e7... Adam's changes to wham and cluster following previous commit #ifdef DEBUG write (iout,*) "The FDIMLESS array before sorting" do i=1,ncon @@ -232,14 +240,12 @@ c write (iout,*) "qfree",qfree sumprob=0.0 do i=1,min0(ncon,maxstr_proc)-1 sumprob=sumprob+dexp(dble(-fdimless(i)+fdimless(1)))/qfree -c#define DEBUG #ifdef DEBUG write (iout,*) 'i=',i,ib,beta_h(ib), & 1.0d0/(1.987d-3*beta_h(ib)),list_conf(i), & totfree(list_conf(i)), & -entfac(list_conf(i)),fdimless(i),sumprob #endif -c#undef DEBUG if (sumprob.gt.prob_limit) goto 122 c if (sumprob.gt.1.00d0) goto 122 nlist=nlist+1 diff --git a/source/cluster/wham/src/readrtns.F b/source/cluster/wham/src/readrtns.F index c40fcbb..2ae3b92 100644 --- a/source/cluster/wham/src/readrtns.F +++ b/source/cluster/wham/src/readrtns.F @@ -36,7 +36,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,'NSTART',nstart,0) call readi(controlcard,'NEND',nend,0) call reada(controlcard,'ECUT',ecut,10.0d0) @@ -45,7 +46,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 @@ -64,6 +66,22 @@ C call readi(controlcard,'CONSTR_DIST',constr_dist,0) write (iout,*) "with_dihed_constr ",with_dihed_constr, & " CONSTR_DIST",constr_dist +<<<<<<< HEAD +======= + + call readi(controlcard,'CONSTR_HOMOL',constr_homology,0) + write (iout,*) "with_homology_constr ",with_dihed_constr, + & " CONSTR_HOMOLOGY",constr_homology + print_homology_restraints= + & index(controlcard,"PRINT_HOMOLOGY_RESTRAINTS").gt.0 + print_contact_map=index(controlcard,"PRINT_CONTACT_MAP").gt.0 + print_homology_models= + & index(controlcard,"PRINT_HOMOLOGY_MODELS").gt.0 + +#ifdef AIX + call flush_(iout) +#else +>>>>>>> 3d6f9e7... Adam's changes to wham and cluster following previous commit call flush(iout) if (min_var) iopt=1 return @@ -240,6 +258,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) @@ -295,6 +326,50 @@ C Convert sequence to numeric code if (nend.gt.nct .or. nend.eq.0) nend=nct write (iout,*) "nstart",nstart," nend",nend nres0=nres +<<<<<<< HEAD +======= + +C Juyong:READ init_vars +C Initialize variables! +C Juyong:READ read_info +C READ fragment information!! +C both routines should be in dfa.F file!! + + if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and. + & wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then +#ifdef DEBUG + write (iout,*) "Calling init_dfa_vars" +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif +#endif + call init_dfa_vars +#ifdef DEBUG + write (iout,*) 'init_dfa_vars finished!' +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif +#endif + call read_dfa_info +#ifdef DEBUG + write (iout,*) 'read_dfa_info finished!' +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif +#endif + endif + + if (constr_homology.gt.0) then + call read_constr_homology(print_homology_restraints) + endif + +>>>>>>> 3d6f9e7... Adam's changes to wham and cluster following previous commit c if (pdbref) then c read(inp,'(a)') pdbfile c write (iout,'(2a)') 'PDB data will be read from file ',pdbfile @@ -362,7 +437,7 @@ c endif enddo enddo endif - call contact(.true.,ncont_ref,icont_ref) + call contact(print_contact_map,ncont_ref,icont_ref) endif c Read distance restraints if (constr_dist.gt.0) then @@ -801,3 +876,411 @@ c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) call flush(iout) return end +<<<<<<< HEAD +======= + +c====------------------------------------------------------------------- + subroutine read_constr_homology(lprn) + + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SETUP' + include 'COMMON.CONTROL' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.GEO' + include 'COMMON.INTERACT' + include 'COMMON.HOMRESTR' +c +c For new homol impl +c + include 'COMMON.VAR' +c include 'include_unres/COMMON.VAR' +c + +c double precision odl_temp,sigma_odl_temp,waga_theta,waga_d, +c & dist_cut +c common /przechowalnia/ odl_temp(maxres,maxres,max_template), +c & sigma_odl_temp(maxres,maxres,max_template) + character*2 kic2 + character*24 model_ki_dist, model_ki_angle + character*500 controlcard + integer ki, i, j, k, l + logical lprn + logical unres_pdb +c +c FP - Nov. 2014 Temporary specifications for new vars +c + double precision rescore_tmp,x12,y12,z12 + double precision, dimension (max_template,maxres) :: rescore + character*24 tpl_k_rescore +c ----------------------------------------------------------------- +c Reading multiple PDB ref structures and calculation of retraints +c not using pre-computed ones stored in files model_ki_{dist,angle} +c FP (Nov., 2014) +c ----------------------------------------------------------------- +c +c +c Alternative: reading from input +#ifdef DEBUG + write (iout,*) "BEGIN READ HOMOLOGY INFO" +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif +#endif + call card_concat(controlcard) + call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0) + call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0) + call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new + call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new + call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma + + call readi(controlcard,"HOMOL_NSET",homol_nset,1) + if (homol_nset.gt.1)then + call readi(controlcard,"ISET",iset,1) + call card_concat(controlcard) + read(controlcard,*) (waga_homology(i),i=1,homol_nset) + else + iset=1 + waga_homology(1)=1.0 + endif +c +#ifdef DEBUG + write(iout,*) "read_constr_homology iset",iset + write(iout,*) "waga_homology(",iset,")",waga_homology(iset) +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif +#endif +cd write (iout,*) "nnt",nnt," nct",nct +cd call flush(iout) + + + lim_odl=0 + lim_dih=0 +c +c New +c + lim_theta=0 + lim_xx=0 +c +c Reading HM global scores (prob not required) +c +c open (4,file="HMscore") +c do k=1,constr_homology +c read (4,*,end=521) hmscore_tmp +c hmscore(k)=hmscore_tmp ! Another transformation can be used +c write(*,*) "Model", k, ":", hmscore(k) +c enddo +c521 continue + +c write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d + + write (iout,*) "CONSTR_HOMOLOGY",constr_homology + do k=1,constr_homology + + read(inp,'(a)') pdbfile +c write (iout,*) "k ",k," pdbfile ",pdbfile +c Next stament causes error upon compilation (?) +c if(me.eq.king.or. .not. out1file) +c write (iout,'(2a)') 'PDB data will be read from file ', +c & pdbfile(:ilen(pdbfile)) + open(ipdbin,file=pdbfile,status='old',err=33) + goto 34 + 33 write (iout,'(a)') 'Error opening PDB file.' + stop + 34 continue +c print *,'Begin reading pdb data' +c +c Files containing res sim or local scores (former containing sigmas) +c + + write(kic2,'(bz,i2.2)') k + + tpl_k_rescore="template"//kic2//".sco" +c tpl_k_sigma_odl="template"//kic2//".sigma_odl" +c tpl_k_sigma_dih="template"//kic2//".sigma_dih" +c tpl_k_sigma_theta="template"//kic2//".sigma_theta" +c tpl_k_sigma_d="template"//kic2//".sigma_d" + + unres_pdb=.false. + call readpdb + do i=1,2*nres + do j=1,3 + crefjlee(j,i)=c(j,i) + enddo + enddo +#ifdef DEBUG + do i=1,nres + write (iout,'(i5,3f8.3,5x,3f8.3)') i,(crefjlee(j,i),j=1,3), + & (crefjlee(j,i+nres),j=1,3) + enddo + write (iout,*) "READ HOMOLOGY INFO" + write (iout,*) "read_constr_homology x: after reading pdb file" + write (iout,*) "waga_homology(",iset,")",waga_homology(iset) + write (iout,*) "waga_dist",waga_dist + write (iout,*) "waga_angle",waga_angle + write (iout,*) "waga_theta",waga_theta + write (iout,*) "waga_d",waga_d + write (iout,*) "dist_cut",dist_cut +#endif +#ifdef AIX + call flush_(iout) +#else + call flush(iout) +#endif + +c +c Distance restraints +c +c ... --> odl(k,ii) +C Copy the coordinates from reference coordinates (?) + do i=1,2*nres + do j=1,3 + c(j,i)=cref(j,i) +c write (iout,*) "c(",j,i,") =",c(j,i) + enddo + enddo +c +c From read_dist_constr (commented out 25/11/2014 <-> res sim) +c +c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore + open (ientin,file=tpl_k_rescore,status='old') + do irec=1,maxdim ! loop for reading res sim + if (irec.eq.1) then + rescore(k,irec)=0.0d0 + goto 1301 + endif + read (ientin,*,end=1401) rescore_tmp +c rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values + rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores +c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec) + 1301 continue + enddo + 1401 continue + close (ientin) +c open (ientin,file=tpl_k_sigma_odl,status='old') +c do irec=1,maxdim ! loop for reading sigma_odl +c read (ientin,*,end=1401) i, j, +c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?) +c sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose? +c & sigma_odl_temp(i+nnt-1,j+nnt-1,k) +c enddo +c 1401 continue +c close (ientin) + if (waga_dist.ne.0.0d0) then + ii=0 + do i = nnt,nct-2 ! right? without parallel. + do j=i+2,nct ! right? +c do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology +c do j=i+2,nres ! ibid +c do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above) +c do j=i+2,nct ! ibid + ii=ii+1 +c write (iout,*) "k",k +c write (iout,*) "i",i," j",j," constr_homology", +c & constr_homology + ires_homo(ii)=i + jres_homo(ii)=j +c +c Attempt to replace dist(i,j) by its definition in ... +c + x12=c(1,i)-c(1,j) + y12=c(2,i)-c(2,j) + z12=c(3,i)-c(3,j) + distal=dsqrt(x12*x12+y12*y12+z12*z12) + odl(k,ii)=distal +c +c odl(k,ii)=dist(i,j) +c write (iout,*) "dist(",i,j,") =",dist(i,j) +c write (iout,*) "distal = ",distal +c write (iout,*) "odl(",k,ii,") =",odl(k,ii) +c write(iout,*) "rescore(",k,i,") =",rescore(k,i), +c & "rescore(",k,j,") =",rescore(k,j) +c +c Calculation of sigma from res sim +c +c if (odl(k,ii).le.6.0d0) then +c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j) +c Other functional forms possible depending on odl(k,ii), eg. +c + if (odl(k,ii).le.dist_cut) then + sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible +c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j) + else +#ifdef OLDSIGMA + sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error + & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) +#else + sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error + & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) +#endif +c Following expr replaced by a positive exp argument +c sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* +c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2) + +c sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)* +c & dexp(-0.5d0*(odl(k,ii)/dist_cut)**2) + endif +c + sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error +c sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii) +c +c sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?) +c & sigma_odl_temp(i,j,k) ! not inverse because of use of res. similarity + enddo +c read (ientin,*) sigma_odl(k,ii) ! 1st variant + enddo +c lim_odl=ii +c if (constr_homology.gt.0) call homology_partition + endif +c +c Theta, dihedral and SC retraints +c + if (waga_angle.gt.0.0d0) then +c open (ientin,file=tpl_k_sigma_dih,status='old') +c do irec=1,maxres-3 ! loop for reading sigma_dih +c read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for? +c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right? +c sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity +c & sigma_dih(k,i+nnt-1) +c enddo +c1402 continue +c close (ientin) + do i = nnt+3,nct ! right? without parallel. +c do i=1,nres ! alternative for bounds acc to readpdb? +c do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology +c do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel. + dih(k,i)=phiref(i) ! right? +c read (ientin,*) sigma_dih(k,i) ! original variant +c write (iout,*) "dih(",k,i,") =",dih(k,i) +c write(iout,*) "rescore(",k,i,") =",rescore(k,i), +c & "rescore(",k,i-1,") =",rescore(k,i-1), +c & "rescore(",k,i-2,") =",rescore(k,i-2), +c & "rescore(",k,i-3,") =",rescore(k,i-3) + + sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+ + & rescore(k,i-2)+rescore(k,i-3) ! right expression ? +c +c write (iout,*) "Raw sigmas for dihedral angle restraints" +c write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i) +c sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* +c rescore(k,i-2)*rescore(k,i-3) ! right expression ? +c Instead of res sim other local measure of b/b str reliability possible + sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i)) +c sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i) + if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right? +c if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file + enddo + endif + + if (waga_theta.gt.0.0d0) then +c open (ientin,file=tpl_k_sigma_theta,status='old') +c do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds? +c read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for? +c sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity +c & sigma_theta(k,i+nnt-1) +c enddo +c1403 continue +c close (ientin) + + do i = nnt+2,nct ! right? without parallel. +c do i = i=1,nres ! alternative for bounds acc to readpdb? +c do i=ithet_start,ithet_end ! with FG parallel. + thetatpl(k,i)=thetaref(i) +c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i) +c write(iout,*) "rescore(",k,i,") =",rescore(k,i), +c & "rescore(",k,i-1,") =",rescore(k,i-1), +c & "rescore(",k,i-2,") =",rescore(k,i-2) +c read (ientin,*) sigma_theta(k,i) ! 1st variant + sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+ + & rescore(k,i-2) ! right expression ? + sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) + +c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* +c rescore(k,i-2) ! right expression ? +c sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i) + if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right? + enddo + endif + + if (waga_d.gt.0.0d0) then +c open (ientin,file=tpl_k_sigma_d,status='old') +c do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds? +c read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for? +c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity +c & sigma_d(k,i+nnt-1) +c enddo +c1404 continue + close (ientin) + + do i = nnt,nct ! right? without parallel. +c do i=2,nres-1 ! alternative for bounds acc to readpdb? +c do i=loc_start,loc_end ! with FG parallel. + if (itype(i).eq.10) goto 1 ! right? + xxtpl(k,i)=xxref(i) + yytpl(k,i)=yyref(i) + zztpl(k,i)=zzref(i) +c write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i) +c write (iout,*) "yytpl(",k,i,") =",yytpl(k,i) +c write (iout,*) "zztpl(",k,i,") =",zztpl(k,i) +c write(iout,*) "rescore(",k,i,") =",rescore(k,i) + sigma_d(k,i)=rescore(k,i) ! right expression ? + sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i)) + +c sigma_d(k,i)=hmscore(k)*rescore(k,i) ! right expression ? +c sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i) +c read (ientin,*) sigma_d(k,i) ! 1st variant + if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right? + 1 continue + enddo + endif + close(ientin) + enddo + if (waga_dist.ne.0.0d0) lim_odl=ii + if (constr_homology.gt.0) call homology_partition + if (constr_homology.gt.0) call init_int_table +cd write (iout,*) "homology_partition: lim_theta= ",lim_theta, +cd & "lim_xx=",lim_xx +c write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end +c write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end +c +c Print restraints +c + if (.not.lprn) return +cd write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d + if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then + write (iout,*) "Distance restraints from templates" + do ii=1,lim_odl + write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii), + & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology) + enddo + write (iout,*) "Dihedral angle restraints from templates" + do i=nnt+3,lim_dih + write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i), + & rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology) + enddo + write (iout,*) "Virtual-bond angle restraints from templates" + do i=nnt+2,lim_theta + write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i), + & rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology) + enddo + write (iout,*) "SC restraints from templates" + do i=nnt,lim_xx + write(iout,'(i5,10(4f8.2,4x))') i, + & (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i), + & 1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology) + enddo + endif +c ----------------------------------------------------------------- + return + end + + +>>>>>>> 3d6f9e7... Adam's changes to wham and cluster following previous commit diff --git a/source/cluster/wham/src/wrtclust.f b/source/cluster/wham/src/wrtclust.f index 97592b7..3915ebc 100644 --- a/source/cluster/wham/src/wrtclust.f +++ b/source/cluster/wham/src/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 diff --git a/source/wham/src/energy_p_new.F b/source/wham/src/energy_p_new.F index 13fe796..2dc4aa7 100644 --- a/source/wham/src/energy_p_new.F +++ b/source/wham/src/energy_p_new.F @@ -4601,6 +4601,12 @@ c 3 = SC...Ca...Ca...SCi cosphi=dcos(j*tauangle(intertyp,i)) sinphi=dsin(j*tauangle(intertyp,i)) esccor=esccor+v1ij*cosphi+v2ij*sinphi +<<<<<<< HEAD +======= +#ifdef DEBUG + esccor_ii=esccor_ii+v1ij*cosphi+v2ij*sinphi +#endif +>>>>>>> 3d6f9e7... Adam's changes to wham and cluster following previous commit gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi) enddo gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci diff --git a/source/wham/src/molread_zs.F b/source/wham/src/molread_zs.F index 6e0727f..534e6df 100644 --- a/source/wham/src/molread_zs.F +++ b/source/wham/src/molread_zs.F @@ -48,6 +48,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 write (iout,*) "Numeric code:" write (iout,'(20i4)') (itype(i),i=1,nres) do i=1,nres-1 -- 1.7.9.5