X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Freadrtns.F;h=3aa5fd6d1bc0cd38cc0af2f5eb5cc3533d687d51;hb=4c8b3525a22636eb28709e3f43c91127f57f8637;hp=056bfd43c91cba638c14e21ae33bc7310187ed72;hpb=c3d632cf84fd647434c093cb3fd69ea605293f75;p=unres.git diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F index 056bfd4..3aa5fd6 100644 --- a/source/unres/src_MD/readrtns.F +++ b/source/unres/src_MD/readrtns.F @@ -8,6 +8,7 @@ include 'COMMON.CONTROL' include 'COMMON.SBRIDGE' include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' logical file_exist C Read force-field parameters except weights call parmread @@ -130,6 +131,7 @@ C Set up the time limit (caution! The time must be input in minutes!) sideadd=(index(controlcard,'SIDEADD').gt.0) energy_dec=(index(controlcard,'ENERGY_DEC').gt.0) outpdb=(index(controlcard,'PDBOUT').gt.0) + outx=(index(controlcard,'XOUT').gt.0) outmol2=(index(controlcard,'MOL2OUT').gt.0) pdbref=(index(controlcard,'PDBREF').gt.0) refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0) @@ -386,6 +388,11 @@ C large = index(controlcard,"LARGE").gt.0 print_compon = index(controlcard,"PRINT_COMPON").gt.0 rattle = index(controlcard,"RATTLE").gt.0 + preminim = index(controlcard,"PREMINIM").gt.0 + if (preminim) then + dccart=(index(controlcard,'CART').gt.0) + call read_minim + endif c if performing umbrella sampling, fragments constrained are read from the fragment file nset=0 if(usampl) then @@ -426,6 +433,10 @@ c if performing umbrella sampling, fragments constrained are read from the frag write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx if (rattle) write (iout,'(a60)') & "Rattle algorithm used to constrain the virtual bonds" + if (preminim .or. iranconf.gt.0) then + write (iout,'(a60)') + & "Initial structure will be energy-minimized" + endif endif reset_fricmat=1000 if (lang.gt.0) then @@ -627,6 +638,7 @@ C common /pizda/ itype_pdb logical seq_comp,fail double precision energia(0:n_ene) + integer ilen external ilen C @@ -922,6 +934,17 @@ C 12/1/95 Added weight for the multi-body term WCORR 34 continue c print *,'Begin reading pdb data' 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 +#endif c print *,'Finished reading pdb data' if(me.eq.king.or..not.out1file) & write (iout,'(a,i3,a,i3)')'nsup=',nsup, @@ -1064,6 +1087,7 @@ C Juyong:READ read_info C READ fragment information!! C both routines should be in dfa.F file!! +#ifdef DFA 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 @@ -1071,6 +1095,7 @@ C both routines should be in dfa.F file!! call read_dfa_info print*, 'read_dfa_info finished!' endif +#endif C if (pdbref) then if(me.eq.king.or..not.out1file) @@ -1165,6 +1190,44 @@ c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup if (constr_homology.gt.0) then call read_constr_homology + if (indpdb.gt.0 .or. pdbref) then + do i=1,2*nres + do j=1,3 + c(j,i)=crefjlee(j,i) + cref(j,i)=crefjlee(j,i) + enddo + enddo + endif +#ifdef DEBUG + write (iout,*) "Array C" + do i=1,nres + write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3), + & (c(j,i+nres),j=1,3) + enddo + write (iout,*) "Array Cref" + do i=1,nres + write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3), + & (cref(j,i+nres),j=1,3) + enddo +#endif + call int_from_cart1(.false.) + call sc_loc_geom(.false.) + do i=1,nres + thetaref(i)=theta(i) + phiref(i)=phi(i) + enddo + 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 else homol_nset=0 endif @@ -1182,6 +1245,8 @@ C initial geometry. if(me.eq.king.or..not.out1file .and.fg_rank.eq.0) & write (iout,'(a)') 'Initial geometry will be read in.' if (read_cart) then + read (inp,*) time,potE,uconst,t_bath, + & nss,(ihpb(j),jhpb(j),j=1,nss), nn, (qfrag(i),i=1,nn) read(inp,'(8f10.5)',end=36,err=36) & ((c(l,k),l=1,3),k=1,nres), & ((c(l,k+nres),l=1,3),k=nnt,nct) @@ -1200,7 +1265,7 @@ C initial geometry. enddo endif enddo - return +c return else call read_angles(inp,*36) endif @@ -1284,6 +1349,8 @@ C Generate distance constraints, if the PDB structure is to be regularized. if (nthread.gt.0) then call read_threadbase endif + write (iout,*) "READRTNS: Calling setup_var" + call flush(iout) call setup_var if (me.eq.king .or. .not. out1file) & call intout @@ -2262,6 +2329,8 @@ c print *,"Processor",myrank," fg_rank",fg_rank & //'.pdb' mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// & liczba(:ilen(liczba))//'.mol2' + cartname=prefix(:lenpre)//'_'//pot(:lenpot)// + & liczba(:ilen(liczba))//'.x' statname=prefix(:lenpre)//'_'//pot(:lenpot)// & liczba(:ilen(liczba))//'.stat' if (lentmp.gt.0) @@ -2280,6 +2349,7 @@ c print *,"Processor",myrank," fg_rank",fg_rank intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int' pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb' mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2' + cartname=prefix(:lenpre)//'_'//pot(:lenpot)//'.x' statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat' if (lentmp.gt.0) & call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) @@ -2440,33 +2510,36 @@ c------------------------------------------------------------------------------- include 'COMMON.IOUNITS' include 'COMMON.MD' include 'COMMON.CONTROL' + integer iset1 read(inp,*) nset,nfrag,npair,nfrag_back if(me.eq.king.or..not.out1file) & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair, & " nfrag_back",nfrag_back - do iset=1,nset - read(inp,*) mset(iset) + do iset1=1,nset + read(inp,*) mset(iset1) do i=1,nfrag - read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset), - & qinfrag(i,iset) + read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1), + & qinfrag(i,iset1) if(me.eq.king.or..not.out1file) - & write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset), - & ifrag(2,i,iset), qinfrag(i,iset) + & write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1), + & ifrag(2,i,iset1), qinfrag(i,iset1) enddo do i=1,npair - read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset), - & qinpair(i,iset) + read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1), + & qinpair(i,iset1) if(me.eq.king.or..not.out1file) - & write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset), - & ipair(2,i,iset), qinpair(i,iset) + & write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1), + & ipair(2,i,iset1), qinpair(i,iset1) enddo do i=1,nfrag_back - read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset), - & wfrag_back(3,i,iset), - & ifrag_back(1,i,iset),ifrag_back(2,i,iset) + read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1), + & wfrag_back(3,i,iset1), + & ifrag_back(1,i,iset1),ifrag_back(2,i,iset1) if(me.eq.king.or..not.out1file) - & write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset), - & wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset) + & write(iout,*) "A",i,wfrag_back(1,i,iset1), + & wfrag_back(2,i,iset1), + & wfrag_back(3,i,iset1),ifrag_back(1,i,iset1), + & ifrag_back(2,i,iset1) enddo enddo return @@ -2651,8 +2724,9 @@ c & sigma_odl_temp(maxres,maxres,max_template) c c FP - Nov. 2014 Temporary specifications for new vars c - double precision rescore_tmp,x12,y12,z12 + double precision rescore_tmp,x12,y12,z12,rescore2_tmp double precision, dimension (max_template,maxres) :: rescore + double precision, dimension (max_template,maxres) :: rescore2 character*24 pdbfile,tpl_k_rescore c ----------------------------------------------------------------- c Reading multiple PDB ref structures and calculation of retraints @@ -2670,24 +2744,24 @@ c Alternative: reading from input call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma call readi(controlcard,"HOMOL_NSET",homol_nset,1) + read2sigma=(index(controlcard,'READ2SIGMA').gt.0) if (homol_nset.gt.1)then call card_concat(controlcard) - read(controlcard,*) (waga_dist1(i),i=1,homol_nset) - call card_concat(controlcard) - read(controlcard,*) (waga_angle1(i),i=1,homol_nset) - write(iout,*) "iset distance_weight angle_weight" - do i=1,homol_nset - write(iout,*) i,waga_dist1(i),waga_angle1(i) - enddo + read(controlcard,*) (waga_homology(i),i=1,homol_nset) + if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then + write(iout,*) "iset homology_weight " + do i=1,homol_nset + write(iout,*) i,waga_homology(i) + enddo + endif iset=mod(kolor,homol_nset)+1 else iset=1 -c call reada(controlcard,"HOMOL_DIST",waga_dist(1),1.0d0) -c call reada(controlcard,"HOMOL_ANGLE",waga_angle(1),1.0d0) + waga_homology(1)=1.0 endif - write (iout,*) "nnt",nnt," nct",nct - call flush(iout) +cd write (iout,*) "nnt",nnt," nct",nct +cd call flush(iout) lim_odl=0 @@ -2698,30 +2772,8 @@ 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 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)) - 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 @@ -2729,23 +2781,7 @@ 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" - call readpdb -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 @@ -2756,24 +2792,54 @@ c write(iout,*) "tpl_k_rescore - ",tpl_k_rescore rescore(k,irec)=0.0d0 goto 1301 endif - read (ientin,*,end=1401) rescore_tmp + if (read2sigma) then + read (ientin,*,end=1401) rescore2_tmp,rescore_tmp, + & idomain(k,irec) + rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) + rescore2(k,irec)=0.5d0*(rescore2_tmp+0.5d0) + else + idomain(k,irec)=1 + 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 + rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores c write(iout,*) "rescore(",k,irec,") =",rescore(k,irec) + endif 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.gt.0.0d0) then + + 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 + + + + if (waga_dist.ne.0.0d0) then ii=0 do i = nnt,nct-2 ! right? without parallel. do j=i+2,nct ! right? @@ -2781,6 +2847,9 @@ c do i = 1,nres ! alternative for bounds as used to set initial values i 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 + + if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0) then + ii=ii+1 c write (iout,*) "k",k c write (iout,*) "i",i," j",j," constr_homology", @@ -2809,31 +2878,44 @@ 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 + 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 & 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 + 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) - 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 + endif 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 @@ -2852,6 +2934,10 @@ c close (ientin) 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. + if (idomain(k,i).eq.0) then + sigma_dih(k,i)=0.0 + cycle + endif dih(k,i)=phiref(i) ! right? c read (ientin,*) sigma_dih(k,i) ! original variant c write (iout,*) "dih(",k,i,") =",dih(k,i) @@ -2888,6 +2974,10 @@ 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. + if (idomain(k,i).eq.0) then + sigma_theta(k,i)=0.0 + cycle + endif thetatpl(k,i)=thetaref(i) c write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i) c write(iout,*) "rescore(",k,i,") =",rescore(k,i), @@ -2913,12 +3003,16 @@ 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 - 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? + if (itype(i).eq.10) cycle + if (idomain(k,i).eq.0 ) then + sigma_d(k,i)=0.0 + cycle + endif xxtpl(k,i)=xxref(i) yytpl(k,i)=yyref(i) zztpl(k,i)=zzref(i) @@ -2933,45 +3027,44 @@ 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.gt.0.0d0) lim_odl=ii + 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 - write (iout,*) "homology_partition: lim_theta= ",lim_theta, - & "lim_xx=",lim_xx +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 - write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d - write (iout,*) "Distance restraints from templates" - do ii=1,lim_odl +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 + 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 + 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 + 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 - + enddo + endif c ----------------------------------------------------------------- return end