X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Freadrtns_CSA.F;fp=source%2Funres%2Fsrc_MD-M%2Freadrtns_CSA.F;h=91d7c0783d9f86e92191779705b1784fb11d60ab;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=6ab2a4e92beb525bc04d6bf01117d80088d497e6;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 6ab2a4e..91d7c07 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -98,6 +98,10 @@ 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) +C this variable with_theta_constr is the variable which allow to read and execute the +C constrains on theta angles WITH_THETA_CONSTR is the keyword + with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0 + write (iout,*) "with_theta_constr ",with_theta_constr call readi(controlcard,'SYM',symetr,1) call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0 @@ -571,7 +575,7 @@ C integer rescode double precision x(maxvar) character*256 pdbfile - character*320 weightcard + character*400 weightcard character*80 weightcard_t,ucase dimension itype_pdb(maxres) common /pizda/ itype_pdb @@ -713,6 +717,14 @@ C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,"V2SS",v2ss,7.61d0) call reada(weightcard,"V3SS",v3ss,13.7d0) call reada(weightcard,"EBR",ebr,-5.50D0) + call reada(weightcard,"ATRISS",atriss,0.301D0) + call reada(weightcard,"BTRISS",btriss,0.021D0) + call reada(weightcard,"CTRISS",ctriss,1.001D0) + call reada(weightcard,"DTRISS",dtriss,1.001D0) + write (iout,*) "ATRISS=", atriss + write (iout,*) "BTRISS=", btriss + write (iout,*) "CTRISS=", ctriss + write (iout,*) "DTRISS=", dtriss dyn_ss=(index(weightcard,'DYN_SS').gt.0) do i=1,maxres dyn_ss_mask(i)=.false. @@ -733,7 +745,11 @@ C 12/1/95 Added weight for the multi-body term WCORR v2ss=v2ss*wstrain/wsc v3ss=v3ss*wstrain/wsc else - ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain + if (wstrain.ne.0.0) then + ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain + else + ss_depth=0.0 + endif endif if(me.eq.king.or..not.out1file) then @@ -755,9 +771,11 @@ C 12/1/95 Added weight for the multi-body term WCORR 33 write (iout,'(a)') 'Error opening PDB file.' stop 34 continue -c print *,'Begin reading pdb data' +c write (iout,*) 'Begin reading pdb data' +c call flush(iout) call readpdb -c print *,'Finished reading pdb data' +c write (iout,*) 'Finished reading pdb data' +c call flush(iout) if(me.eq.king.or..not.out1file) & write (iout,'(a,i3,a,i3)')'nsup=',nsup, & ' nstart_sup=',nstart_sup @@ -847,27 +865,78 @@ C 8/13/98 Set limits to generating the dihedral angles enddo read (inp,*) ndih_constr if (ndih_constr.gt.0) then - read (inp,*) ftors - read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr) +C read (inp,*) ftors + read (inp,*) (idih_constr(i),phi0(i),drange(i),ftors(i), + & i=1,ndih_constr) if(me.eq.king.or..not.out1file)then write (iout,*) & 'There are',ndih_constr,' constraints on phi angles.' do i=1,ndih_constr - write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i) + write (iout,'(i5,3f8.3)') idih_constr(i),phi0(i),drange(i), + & ftors(i) enddo endif do i=1,ndih_constr phi0(i)=deg2rad*phi0(i) drange(i)=deg2rad*drange(i) enddo - if(me.eq.king.or..not.out1file) - & write (iout,*) 'FTORS',ftors +C if(me.eq.king.or..not.out1file) +C & write (iout,*) 'FTORS',ftors do i=1,ndih_constr ii = idih_constr(i) phibound(1,ii) = phi0(i)-drange(i) phibound(2,ii) = phi0(i)+drange(i) enddo endif +C first setting the theta boundaries to 0 to pi +C this mean that there is no energy penalty for any angle occuring this can be applied +C for generate random conformation but is not implemented in this way +C do i=1,nres +C thetabound(1,i)=0 +C thetabound(2,i)=pi +C enddo +C begin reading theta constrains this is quartic constrains allowing to +C have smooth second derivative + if (with_theta_constr) then +C with_theta_constr is keyword allowing for occurance of theta constrains + read (inp,*) ntheta_constr +C ntheta_constr is the number of theta constrains + if (ntheta_constr.gt.0) then +C read (inp,*) ftors + read (inp,*) (itheta_constr(i),theta_constr0(i), + & theta_drange(i),for_thet_constr(i), + & i=1,ntheta_constr) +C the above code reads from 1 to ntheta_constr +C itheta_constr(i) residue i for which is theta_constr +C theta_constr0 the global minimum value +C theta_drange is range for which there is no energy penalty +C for_thet_constr is the force constant for quartic energy penalty +C E=k*x**4 + if(me.eq.king.or..not.out1file)then + write (iout,*) + & 'There are',ntheta_constr,' constraints on phi angles.' + do i=1,ntheta_constr + write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i), + & theta_drange(i), + & for_thet_constr(i) + enddo + endif + do i=1,ntheta_constr + theta_constr0(i)=deg2rad*theta_constr0(i) + theta_drange(i)=deg2rad*theta_drange(i) + enddo +C if(me.eq.king.or..not.out1file) +C & write (iout,*) 'FTORS',ftors +C do i=1,ntheta_constr +C ii = itheta_constr(i) +C thetabound(1,ii) = phi0(i)-drange(i) +C thetabound(2,ii) = phi0(i)+drange(i) +C enddo + endif ! ntheta_constr.gt.0 + endif! with_theta_constr +C +C with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0 +C write (iout,*) "with_dihed_constr ",with_dihed_constr nnt=1 #ifdef MPI if (me.eq.king) then @@ -954,7 +1023,9 @@ czscore call geom_to_var(nvar,coord_exp_zs(1,1)) enddo call contact(.true.,ncont_ref,icont_ref,co) endif -c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup + endif + print *, "A TU" + write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup call flush(iout) if (constr_dist.gt.0) call read_dist_constr write (iout,*) "After read_dist_constr nhpb",nhpb @@ -975,7 +1046,7 @@ c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup & restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i) enddo endif - endif +C endif if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 & .and. modecalc.ne.8 .and. modecalc.ne.9 .and. & modecalc.ne.10) then @@ -1147,6 +1218,7 @@ cd write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar) & write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)') & 'Processor',myrank,': end reading molecular data.' #endif + print *,"A TU?" return end c-------------------------------------------------------------------------- @@ -2303,7 +2375,8 @@ c------------------------------------------------------------------------------- integer ifrag_(2,100),ipair_(2,100) double precision wfrag_(100),wpair_(100) character*500 controlcard -c write (iout,*) "Calling read_dist_constr" + print *, "WCHODZE" + write (iout,*) "Calling read_dist_constr" c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup c call flush(iout) call card_concat(controlcard) @@ -2397,11 +2470,30 @@ c write (iout,*) "j",j," k",k enddo endif enddo + print *,ndist_ do i=1,ndist_ - read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) + if (constr_dist.eq.11) then + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), + & ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1) + fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) + else +C print *,"in else" + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), + & ibecarb(i),forcon(nhpb+1) + endif if (forcon(nhpb+1).gt.0.0d0) then nhpb=nhpb+1 - dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + if (ibecarb(i).gt.0) then + ihpb(i)=ihpb(i)+nres + jhpb(i)=jhpb(i)+nres + endif + if (dhpb(nhpb).eq.0.0d0) + & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) + endif +C read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1) +C if (forcon(nhpb+1).gt.0.0d0) then +C nhpb=nhpb+1 +C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) #ifdef MPI if (.not.out1file .or. me.eq.king) & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", @@ -2410,7 +2502,7 @@ c write (iout,*) "j",j," k",k write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) #endif - endif + enddo call flush(iout) return