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=ebb59da9d7e43553aac1a1c27d2e5858f10c35aa;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=889e5a59e722fff128d0cebf28d150eada2a4906;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 889e5a5..ebb59da 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -81,6 +81,7 @@ C include 'COMMON.INTERACT' include 'COMMON.SETUP' include 'COMMON.SPLITELE' + include 'COMMON.SHIELD' COMMON /MACHSW/ KDIAG,ICORFL,IXDR character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/ character*80 ucase @@ -98,6 +99,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 @@ -141,6 +146,15 @@ C Set up the time limit (caution! The time must be input in minutes!) selfguide=(index(controlcard,'SELFGUIDE')) print *,'AFMlog',AFMlog,selfguide,"KUPA" call readi(controlcard,'IPRINT',iprint,0) +C SHIELD keyword sets if the shielding effect of side-chains is used +C 0 denots no shielding is used all peptide are equally despite the +C solvent accesible area +C 1 the newly introduced function +C 2 reseved for further possible developement + call readi(controlcard,'SHIELD',shield_mode,0) +C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then + write(iout,*) "shield_mode",shield_mode +C endif call readi(controlcard,'MAXGEN',maxgen,10000) call readi(controlcard,'MAXOVERLAP',maxoverlap,1000) call readi(controlcard,"KDIAG",kdiag,0) @@ -244,8 +258,24 @@ C endif write(iout,*) "bordlipbot=",bordlipbot write(iout,*) "bufliptop=",bufliptop write(iout,*) "buflipbot=",buflipbot - - + write (iout,*) "SHIELD MODE",shield_mode + if (shield_mode.gt.0) then + pi=3.141592d0 +C VSolvSphere the volume of solving sphere +C print *,pi,"pi" +C rpp(1,1) is the energy r0 for peptide group contact and will be used for it +C there will be no distinction between proline peptide group and normal peptide +C group in case of shielding parameters + VSolvSphere=4.0/3.0*pi*rpp(1,1)**3 + VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3 + write (iout,*) VSolvSphere,VSolvSphere_div +C long axis of side chain + do i=1,ntyp + long_r_sidechain(i)=vbldsc0(1,i) + short_r_sidechain(i)=sigma0(i) + enddo + buff_shield=1.0d0 + endif if (me.eq.king .or. .not.out1file ) & write (iout,*) "DISTCHAINMAX",distchainmax @@ -571,7 +601,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 +743,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 +771,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 +797,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 +891,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 +1049,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 +1072,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 @@ -1035,6 +1132,7 @@ C initial geometry. omeg(i)=-120d0*deg2rad if (itype(i).le.0) omeg(i)=-omeg(i) enddo + call chainbuild_extconf else if(me.eq.king.or..not.out1file) & write (iout,'(a)') 'Random-generated initial geometry.' @@ -1147,6 +1245,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 +2402,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 +2497,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 +2529,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