X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Freadrtns.F;h=887e89fd80401ac26255e6f14e70fc7918c3a85a;hb=98ab6f62b288476b57cf515d4d04ddafda983979;hp=73ffe83c721c92f3590f0983cf572d4e3959f5b2;hpb=fb112c2c0d337f5433d28f2c66ecf420102e63f0;p=unres.git diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F index 73ffe83..887e89f 100644 --- a/source/unres/src_MD/readrtns.F +++ b/source/unres/src_MD/readrtns.F @@ -2721,6 +2721,8 @@ c & sigma_odl_temp(maxres,maxres,max_template) character*500 controlcard integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp logical lprn /.true./ + integer ilen + external ilen c c FP - Nov. 2014 Temporary specifications for new vars c @@ -2790,6 +2792,20 @@ c do k=1,constr_homology + read(inp,'(a)') 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)) + write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file', + & pdbfile(:ilen(pdbfile)) + open(ipdbin,file=pdbfile,status='old',err=33) + goto 34 + 33 write (iout,'(a,5x,a)') 'Error opening PDB file', + & pdbfile(:ilen(pdbfile)) + stop + 34 continue +c print *,'Begin reading pdb data' c c Files containing res sim or local scores (former containing sigmas) c @@ -2798,6 +2814,23 @@ c tpl_k_rescore="template"//kic2//".sco" + unres_pdb=.false. + if (read2sigma) then + call readpdb_template(k) + else + call readpdb + 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 @@ -2822,45 +2855,11 @@ c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec) endif enddo 1401 continue - close (ientin) - - read(inp,'(a)') 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' - - unres_pdb=.false. - if (read2sigma) then - call readpdb_template(k) - else - call readpdb - endif -c -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 - - - + 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 + 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 @@ -2881,49 +2880,34 @@ c 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 (read2sigma) then - sigma_odl(k,ii)=0 - do ik=i,j - sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik) - enddo - sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1) - if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = - & 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) - else + if (read2sigma) then + sigma_odl(k,ii)=0 + do ik=i,j + sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik) + enddo + sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1) + if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) = + & 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) + 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))* ! 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 + 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 - endif - 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) @@ -2951,10 +2935,7 @@ 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. + do i = nnt+3,nct if (idomain(k,i).eq.0) then sigma_dih(k,i)=0.0 cycle @@ -2967,9 +2948,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))/4.0 -c + 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 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)* @@ -3004,8 +2985,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)/3.0 + 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)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i)) c sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)* @@ -3023,7 +3005,6 @@ c sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use c & sigma_d(k,i+nnt-1) c enddo c1404 continue - do i = nnt,nct ! right? without parallel. c do i=2,nres-1 ! alternative for bounds acc to readpdb? @@ -3046,6 +3027,7 @@ c write(iout,*) "rescore(",k,i,") =",rescore(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? enddo lim_xx=nct-nnt+1 endif