From 3c5e1a8829ff84c0725642015445f4debd4ecc0a Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Tue, 10 May 2016 20:25:23 +0200 Subject: [PATCH] updates from src_MD to src_MD-M for unres --- source/unres/src_MD-M/COMMON.CHAIN | 4 +- source/unres/src_MD-M/COMMON.CONTROL | 10 ++-- source/unres/src_MD-M/DIMENSIONS | 2 +- source/unres/src_MD-M/MD_A-MTS.F | 24 +++++++++ source/unres/src_MD-M/energy_p_new_barrier.F | 13 ++++- source/unres/src_MD-M/readpdb.F | 1 + source/unres/src_MD-M/readrtns_CSA.F | 74 +++++++++++++------------- 7 files changed, 82 insertions(+), 46 deletions(-) diff --git a/source/unres/src_MD-M/COMMON.CHAIN b/source/unres/src_MD-M/COMMON.CHAIN index 84ebf1b..dd53a8e 100644 --- a/source/unres/src_MD-M/COMMON.CHAIN +++ b/source/unres/src_MD-M/COMMON.CHAIN @@ -3,7 +3,7 @@ & afmbeg double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm,t,r, & prod,rt,dc_work,cref,crefjlee,chain_rep,dc_norm2,velAFMconst, - & totTafm + & totTafm,chomo common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2), & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2), & dc_norm2(3,0:maxres2), @@ -24,4 +24,4 @@ common /afm/ forceAFMconst, distafminit,afmend,afmbeg, & velAFMconst, & totTafm - + common /chomo_models/ chomo(3,maxres2+2,max_template) diff --git a/source/unres/src_MD-M/COMMON.CONTROL b/source/unres/src_MD-M/COMMON.CONTROL index 0a868ce..35fb03d 100644 --- a/source/unres/src_MD-M/COMMON.CONTROL +++ b/source/unres/src_MD-M/COMMON.CONTROL @@ -2,19 +2,21 @@ & inprint,i2ndstr,mucadyn,constr_dist,symetr,AFMlog,selfguide, & constr_homology,homol_nset real*8 waga_homology - real*8 waga_dist, waga_angle, waga_theta, waga_d, dist_cut + real*8 waga_dist, waga_angle, waga_theta, waga_d, dist_cut, + & dist2_cut logical minim,refstr,pdbref,outpdb,outmol2,overlapsc,energy_dec, & sideadd,lsecondary,read_cart,unres_pdb, & vdisulf,searchsc,lmuca,dccart,extconf,out1file, - & gnorm_check,gradout,split_ene,read2sigma + & gnorm_check,gradout,split_ene,read2sigma, + & start_from_model common /cntrl/ modecalc,iscode,indpdb,indback,indphi,iranconf, & icheckgrad,minim,i2ndstr,refstr,pdbref,outpdb,outmol2,iprint, & overlapsc,energy_dec,sideadd,lsecondary,read_cart,unres_pdb & ,vdisulf,searchsc,lmuca,dccart,mucadyn,extconf,out1file, & constr_dist,gnorm_check,gradout,split_ene,symetr,AFMlog, & selfguide, - & constr_homology,homol_nset,read2sigma + & constr_homology,homol_nset,read2sigma,start_from_model common /homol/ waga_homology(maxprocs/20), - & waga_dist, waga_angle, waga_theta, waga_d, dist_cut + & waga_dist, waga_angle, waga_theta, waga_d, dist_cut,dist2_cut C... minim = .true. means DO minimization. C... energy_dec = .true. means print energy decomposition matrix diff --git a/source/unres/src_MD-M/DIMENSIONS b/source/unres/src_MD-M/DIMENSIONS index 22c4844..70ce031 100644 --- a/source/unres/src_MD-M/DIMENSIONS +++ b/source/unres/src_MD-M/DIMENSIONS @@ -16,7 +16,7 @@ C Max. number of coarse-grain processors parameter (max_cg_procs=maxprocs) C Max. number of AA residues integer maxres - parameter (maxres=600) + parameter (maxres=800) C Appr. max. number of interaction sites integer maxres2,maxres6,mmaxres2 parameter (maxres2=2*maxres,maxres6=6*maxres) diff --git a/source/unres/src_MD-M/MD_A-MTS.F b/source/unres/src_MD-M/MD_A-MTS.F index 7a9e307..deef97c 100644 --- a/source/unres/src_MD-M/MD_A-MTS.F +++ b/source/unres/src_MD-M/MD_A-MTS.F @@ -1609,6 +1609,30 @@ c Removing the velocity of the center of mass write (iout,*) (vcm(j),j=1,3) endif if (.not.rest) then + if (start_from_model) then + i_model=iran_num(1,constr_homology) + write (iout,*) 'starting from model ',i_model + do i=1,2*nres + do j=1,3 + c(j,i)=chomo(j,i,i_model) + enddo + enddo + + call int_from_cart(.true.,.false.) + call sc_loc_geom(.false.) + do i=1,nres-1 + do j=1,3 + dc(j,i)=c(j,i+1)-c(j,i) + dc_norm(j,i)=dc(j,i)*vbld_inv(i+1) + enddo + enddo + do i=2,nres-1 + do j=1,3 + dc(j,i+nres)=c(j,i+nres)-c(j,i) + dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres) + enddo + enddo + endif call chainbuild write (iout,*) "PREMINIM ",preminim if(iranconf.ne.0 .or. preminim) then diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 79cd1b9..b3c4c22 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -6896,7 +6896,18 @@ c endif enddo - min_odl=minval(distancek) +c min_odl=minval(distancek) + do kk=1,constr_homology + if(l_homo(kk,ii)) then + min_odl=distancek(kk) + exit + endif + enddo + do kk=1,constr_homology + if(l_homo(kk,ii) .and. distancek(kk).lt.min_odl) + & min_odl=distancek(kk) + enddo + c write (iout,* )"min_odl",min_odl #ifdef DEBUG write (iout,*) "ij dij",i,j,dij diff --git a/source/unres/src_MD-M/readpdb.F b/source/unres/src_MD-M/readpdb.F index a504ff7..de172bd 100644 --- a/source/unres/src_MD-M/readpdb.F +++ b/source/unres/src_MD-M/readpdb.F @@ -1006,6 +1006,7 @@ c write (iout,*) "spraw lancuchy",(c(j,i),j=1,3) endif endif do j=1,3 + chomo(j,i,k)=c(j,i) cref(j,i,cou)=c(j,i) cref(j,i+nres,cou)=c(j,i+nres) if (i.le.nres) then diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 7be445d..f944233 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -2580,9 +2580,19 @@ c Alternative: reading from input 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 reada(controlcard,'DIST2_CUT',dist2_cut,9999.0d0) call readi(controlcard,"HOMOL_NSET",homol_nset,1) read2sigma=(index(controlcard,'READ2SIGMA').gt.0) + start_from_model=(index(controlcard,'START_FROM_MODELS').gt.0) + if(.not.read2sigma.and.start_from_model) then + write(iout,*) 'START_FROM_MODELS works only with READ2SIGMA' + start_from_model=.false. + endif + if(start_from_model) write(iout,*) 'START_FROM_MODELS is ON' + if(start_from_model .and. rest) then + write(iout,*) 'START_FROM_MODELS is OFF' + write(iout,*) 'remove restart keyword from input' + endif if (homol_nset.gt.1)then call card_concat(controlcard) read(controlcard,*) (waga_homology(i),i=1,homol_nset) @@ -2679,8 +2689,8 @@ c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore & idomain_tmp i_tmp=i_tmp+nnt-1 idomain(k,i_tmp)=idomain_tmp - rescore(k,i_tmp)=0.5d0*(rescore_tmp+0.5d0) - rescore2(k,i_tmp)=0.5d0*(rescore2_tmp+0.5d0) + rescore(k,i_tmp)=rescore_tmp + rescore2(k,i_tmp)=rescore2_tmp else idomain(k,irec)=1 read (ientin,*,end=1401) rescore_tmp @@ -2697,8 +2707,15 @@ c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec) do i = nnt,nct-2 do j=i+2,nct - if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0) then - + 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) +c write (iout,*) k,i,j,distal,dist2_cut + + if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0 + & .and. distal.le.dist2_cut ) then + ii=ii+1 ii_in_use(ii)=1 l_homo(k,ii)=.true. @@ -2708,13 +2725,6 @@ 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 if (read2sigma) then sigma_odl(k,ii)=0 @@ -2726,30 +2736,18 @@ c & sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) else 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) + sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) else #ifdef OLDSIGMA - sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error + sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2) #else - sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error + sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* & dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0) #endif endif 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) -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 + sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) else ii=ii+1 l_homo(k,ii)=.false. @@ -2784,9 +2782,9 @@ 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) - if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0 + sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+ + & rescore(k,i-2)+rescore(k,i-3))/4.0 +c if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0 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)* @@ -2821,9 +2819,9 @@ 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) - if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0 + sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+ + & rescore(k,i-2))/3.0 +c if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0 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)* @@ -2908,24 +2906,24 @@ 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(2f8.2,1x,l1,4x))') + write(iout,'(3i5,100(2f8.2,1x,l1,4x))') & ii,ires_homo(ii),jres_homo(ii), & (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(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), + write (iout,'(i5,100(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), + write (iout,'(i5,100(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, + write(iout,'(i5,100(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 -- 1.7.9.5