X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Freadrtns.F;h=18e2c2ed4da1496deb1e4b9ecbf709618251c801;hb=a4e5a256e77beab0f13a18e273b46707d587060a;hp=88b7242bd799dafb482492c658dc8a5fea5d8170;hpb=2ff50af0c3c062eea3ed8938b9d9abf88d040b85;p=unres.git diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F index 88b7242..18e2c2e 100644 --- a/source/unres/src_MD/readrtns.F +++ b/source/unres/src_MD/readrtns.F @@ -96,6 +96,7 @@ c print *,"Processor",me," fg_rank",fg_rank," out1file",out1file C Set up the time limit (caution! The time must be input in minutes!) read_cart=index(controlcard,'READ_CART').gt.0 call readi(controlcard,'CONSTR_DIST',constr_dist,0) + call readi(controlcard,'CONSTR_HOMOL',constr_homology,0) call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes @@ -749,6 +750,10 @@ C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,'WTORD',wtor_d,1.0D0) call reada(weightcard,'WANG',wang,1.0D0) call reada(weightcard,'WSCLOC',wscloc,1.0D0) + call reada(weightcard,'WDFAD',wdfa_dist,0.0d0) + call reada(weightcard,'WDFAT',wdfa_tor,0.0d0) + call reada(weightcard,'WDFAN',wdfa_nei,0.0d0) + call reada(weightcard,'WDFAB',wdfa_beta,0.0d0) call reada(weightcard,'SCAL14',scal14,0.4D0) call reada(weightcard,'SCALSCP',scalscp,1.0d0) call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) @@ -778,11 +783,16 @@ C 12/1/95 Added weight for the multi-body term WCORR weights(18)=scal14 weights(21)=wsccor endif + weights(25)=wdfa_dist + weights(26)=wdfa_tor + weights(27)=wdfa_nei + weights(28)=wdfa_beta if(me.eq.king.or..not.out1file) & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3, - & wturn4,wturn6 + & wturn4,wturn6, + & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta 10 format (/'Energy-term weights (unscaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ @@ -801,7 +811,11 @@ C 12/1/95 Added weight for the multi-body term WCORR & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ - & 'WTURN6= ',f10.6,' (turns, 6th order)') + & 'WTURN6= ',f10.6,' (turns, 6th order)'/ + & 'WDFA_D= ',f10.6,' (DFA, distance)' / + & 'WDFA_T= ',f10.6,' (DFA, torsional)' / + & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' / + & 'WDFA_B= ',f10.6,' (DFA, beta formation)') if(me.eq.king.or..not.out1file)then if (wcorr4.gt.0.0d0) then write (iout,'(/2a/)') 'Local-electrostatic type correlation ', @@ -829,7 +843,8 @@ C 12/1/95 Added weight for the multi-body term WCORR if(me.eq.king.or..not.out1file) & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor, & wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3, - & wturn4,wturn6 + & wturn4,wturn6, + & wdfa_dist,wdfa_tor,wdfa_nei,wdfa_beta 22 format (/'Energy-term weights (scaled):'// & 'WSCC= ',f10.6,' (SC-SC)'/ & 'WSCP= ',f10.6,' (SC-p)'/ @@ -848,7 +863,11 @@ C 12/1/95 Added weight for the multi-body term WCORR & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/ & 'WTURN3= ',f10.6,' (turns, 3rd order)'/ & 'WTURN4= ',f10.6,' (turns, 4th order)'/ - & 'WTURN6= ',f10.6,' (turns, 6th order)') + & 'WTURN6= ',f10.6,' (turns, 6th order)'/ + & 'WDFA_D= ',f10.6,' (DFA, distance)' / + & 'WDFA_T= ',f10.6,' (DFA, torsional)' / + & 'WDFA_N= ',f10.6,' (DFA, number of neighbor)' / + & 'WDFA_B= ',f10.6,' (DFA, beta formation)') if(me.eq.king.or..not.out1file) & write (iout,*) "Reference temperature for weights calculation:", & temp0 @@ -940,8 +959,8 @@ C 10/03/12 Adam: Recalculate coordinates with new side chain positions call chainbuild endif C Following 2 lines for diagnostics; comment out if not needed - write (iout,*) "After sideadd" - call intout +c write (iout,*) "After sideadd" +c call intout endif if (indpdb.eq.0) then C Read sequence if not taken from the pdb file. @@ -1039,6 +1058,21 @@ C 8/13/98 Set limits to generating the dihedral angles cd print *,'NNT=',NNT,' NCT=',NCT if (itype(1).eq.21) nnt=2 if (itype(nres).eq.21) nct=nct-1 + +C Bartek: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 + call init_dfa_vars + print*, 'init_dfa_vars finished!' + call read_dfa_info + print*, 'read_dfa_info finished!' + endif +C if (pdbref) then if(me.eq.king.or..not.out1file) & write (iout,'(a,i3)') 'nsup=',nsup @@ -1127,8 +1161,17 @@ czscore call geom_to_var(nvar,coord_exp_zs(1,1)) c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup if (constr_dist.gt.0) then call read_dist_constr - call hpb_partition endif + + + if (constr_homology.gt.0) then + call read_constr_homology + else + homol_nset=0 + endif + + + if (nhpb.gt.0) call hpb_partition c write (iout,*) "After read_dist_constr nhpb",nhpb c call flush(iout) if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 @@ -1251,18 +1294,6 @@ C Generate distance constraints, if the PDB structure is to be regularized. write (iout,'(20i4)') (iss(i),i=1,ns) if (dyn_ss) then write(iout,*)"Running with dynamic disulfide-bond formation" - do i=nss+1,nhpb - ihpb(i-nss)=ihpb(i) - jhpb(i-nss)=jhpb(i) - forcon(i-nss)=forcon(i) - dhpb(i-nss)=dhpb(i) - enddo - nhpb=nhpb-nss - nss=0 - call hpb_partition - do i=1,ns - dyn_ss_mask(iss(i))=.true. - enddo else write (iout,'(/a/)') 'Pre-formed links are:' do i=1,nss @@ -1270,14 +1301,27 @@ C Generate distance constraints, if the PDB structure is to be regularized. i2=jhpb(i)-nres it1=itype(i1) it2=itype(i2) - if (me.eq.king.or..not.out1file) - & write (iout,'(2a,i3,3a,i3,a,3f10.3)') + write (iout,'(2a,i3,3a,i3,a,3f10.3)') & restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i), & ebr,forcon(i) enddo write (iout,'(a)') endif endif + if (ns.gt.0.and.dyn_ss) then + do i=nss+1,nhpb + ihpb(i-nss)=ihpb(i) + jhpb(i-nss)=jhpb(i) + forcon(i-nss)=forcon(i) + dhpb(i-nss)=dhpb(i) + enddo + nhpb=nhpb-nss + nss=0 + call hpb_partition + do i=1,ns + dyn_ss_mask(iss(i))=.true. + enddo + endif if (i2ndstr.gt.0) call secstrp2dihc c call geom_to_var(nvar,x) c call etotal(energia(0)) @@ -1340,10 +1384,12 @@ C Check whether the specified bridging residues are cystines. do i=1,ns if (itype(iss(i)).ne.1) then if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)') - & 'Do you REALLY think that the residue ',restyp(iss(i)),i, + & 'Do you REALLY think that the residue ', + & restyp(itype(iss(i))),i, & ' can form a disulfide bridge?!!!' write (*,'(2a,i3,a)') - & 'Do you REALLY think that the residue ',restyp(iss(i)),i, + & 'Do you REALLY think that the residue ', + & restyp(itype(iss(i))),i, & ' can form a disulfide bridge?!!!' #ifdef MPI call MPI_Finalize(MPI_COMM_WORLD,ierror) @@ -1354,7 +1400,8 @@ C Check whether the specified bridging residues are cystines. C Read preformed bridges. if (ns.gt.0) then read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss) - write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss) + if(fg_rank.eq.0) + & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss) if (nss.gt.0) then nhpb=nss C Check if the residues involved in bridges are in the specified list of @@ -2367,6 +2414,7 @@ c------------------------------------------------------------------------------- include 'COMMON.CHAIN' include 'COMMON.IOUNITS' include 'COMMON.MD' + include 'COMMON.CONTROL' open(irest2,file=rest2name,status='unknown') read(irest2,*) totT,EK,potE,totE,t_bath do i=1,2*nres @@ -2375,7 +2423,7 @@ c------------------------------------------------------------------------------- do i=1,2*nres read(irest2,'(3e15.5)') (dc(j,i),j=1,3) enddo - if(usampl) then + if(usampl.or.homol_nset.gt.1) then read (irest2,*) iset endif close(irest2) @@ -2482,7 +2530,7 @@ c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) if (wfrag_(i).gt.0.0d0) then do j=ifrag_(1,i),ifrag_(2,i)-1 do k=j+1,ifrag_(2,i) - write (iout,*) "j",j," k",k +c write (iout,*) "j",j," k",k ddjk=dist(j,k) if (constr_dist.eq.1) then nhpb=nhpb+1 @@ -2572,6 +2620,135 @@ c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) return end c------------------------------------------------------------------------------- + + subroutine read_constr_homology + + include 'DIMENSIONS' +#ifdef MPI + include 'mpif.h' +#endif + include 'COMMON.SETUP' + include 'COMMON.CONTROL' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.MD' + include 'COMMON.GEO' + include 'COMMON.INTERACT' + double precision odl_temp,sigma_odl_temp + common /przechowalnia/ odl_temp(maxres,maxres,max_template), + & sigma_odl_temp(maxres,maxres,max_template) + character*2 kic2 + character*24 model_ki_dist, model_ki_angle + character*500 controlcard + character*3200 controlcard1 + integer ki, i, j, k, l + logical lprn /.true./ + + call card_concat(controlcard) + call readi(controlcard,"HOMOL_NSET",homol_nset,1) + if (homol_nset.gt.1)then + call card_concat(controlcard) + read(controlcard,*) (waga_dist(i),i=1,homol_nset) + call card_concat(controlcard) + read(controlcard,*) (waga_angle(i),i=1,homol_nset) + write(iout,*) "iset distance_weight angle_weight" + do i=1,homol_nset + write(iout,*) i,waga_dist(i),waga_angle(i) + enddo + else + iset=1 + call reada(controlcard,"HOMOL_DIST",waga_dist(1),1.0d0) + call reada(controlcard,"HOMOL_ANGLE",waga_angle(1),1.0d0) + endif + + write (iout,*) "nnt",nnt," nct",nct + call flush(iout) + lim_odl=0 + lim_dih=0 + do i=1,nres + do j=i+2,nres + do ki=1,constr_homology + sigma_odl_temp(i,j,ki)=0.0d0 + odl_temp(i,j,ki)=0.0d0 + enddo + enddo + enddo + do i=1,nres-3 + do ki=1,constr_homology + dih(ki,i)=0.0d0 + sigma_dih(ki,i)=0.0d0 + enddo + enddo + do ki=1,constr_homology + write(kic2,'(i2)') ki + if (ki.le.9) kic2="0"//kic2(2:2) + + model_ki_dist="model"//kic2//".dist" + model_ki_angle="model"//kic2//".angle" + open (ientin,file=model_ki_dist,status='old') + do irec=1,maxdim !petla do czytania wiezow na odleglosc + read (ientin,*,end=1401) i, j, odl_temp(i+nnt-1,j+nnt-1,ki), + & sigma_odl_temp(i+nnt-1,j+nnt-1,ki) + odl_temp(j+nnt-1,i+nnt-1,ki)=odl_temp(i+nnt-1,j+nnt-1,ki) + sigma_odl_temp(j+nnt-1,i+nnt-1,ki)= + & sigma_odl_temp(i+nnt-1,j+nnt-1,ki) + enddo + 1401 continue + close (ientin) + open (ientin,file=model_ki_angle,status='old') + do irec=1,maxres-3 !petla do czytania wiezow na katy torsyjne + read (ientin,*,end=1402) i, j, k,l,dih(ki,i+nnt-1), + & sigma_dih(ki,i+nnt-1) + if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 + sigma_dih(ki,i+nnt-1)=1.0d0/sigma_dih(ki,i+nnt-1)**2 + enddo + 1402 continue + close (ientin) + enddo + ii=0 + write (iout,*) "nnt",nnt," nct",nct + do i=nnt,nct-2 + do j=i+2,nct + ki=1 +c write (iout,*) "i",i," j",j," constr_homology",constr_homology + do while (ki.le.constr_homology .and. + & sigma_odl_temp(i,j,ki).le.0.0d0) +c write (iout,*) "ki",ki," sigma_odl",sigma_odl_temp(i,j,ki) + ki=ki+1 + enddo +c write (iout,*) "ki",ki + if (ki.gt.constr_homology) cycle + ii=ii+1 + ires_homo(ii)=i + jres_homo(ii)=j + do ki=1,constr_homology + odl(ki,ii)=odl_temp(i,j,ki) + sigma_odl(ki,ii)=1.0d0/sigma_odl_temp(i,j,ki)**2 + enddo + enddo + enddo + lim_odl=ii + if (constr_homology.gt.0) call homology_partition +c Print restraints + if (.not.lprn) return + write (iout,*) "Distance restraints from templates" + do ii=1,lim_odl + write(iout,'(3i5,10(2f8.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,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 +c write(iout,*) "TEST CZYTANIA1",odl(1,2,1),odl(1,3,1),odl(1,4,1) +c write(iout,*) "TEST CZYTANIA2",dih(1,1),dih(2,1),dih(3,1) + + + return + end +c---------------------------------------------------------------------- + #ifdef WINIFL subroutine flush(iu) return @@ -2583,6 +2760,7 @@ c------------------------------------------------------------------------------- return end #endif + c------------------------------------------------------------------------------ subroutine copy_to_tmp(source) include "DIMENSIONS"