X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-M%2Freadrtns.F;h=1bf78e886146532d8a9d471b40334877f0ec309f;hb=fd10e81284102e5a8466e05c87f0202979a2f00b;hp=36c13b17068792d185c630bf54a39c5de53d2a4b;hpb=478a9d9a1c99eb3f4bc4ca676ff3162bdd01d633;p=unres.git diff --git a/source/wham/src-M/readrtns.F b/source/wham/src-M/readrtns.F index 36c13b1..1bf78e8 100644 --- a/source/wham/src-M/readrtns.F +++ b/source/wham/src-M/readrtns.F @@ -17,6 +17,9 @@ include "COMMON.FREE" include "COMMON.CONTROL" include "COMMON.ENERGIES" + include "COMMON.SPLITELE" + include "COMMON.SBRIDGE" + include "COMMON.SHIELD" character*800 controlcard integer i,j,k,ii,n_ene_found integer ind,itype1,itype2,itypf,itypsc,itypp @@ -25,7 +28,7 @@ character*16 ucase character*16 key external ucase - + double precision pi call card_concat(controlcard,.true.) call readi(controlcard,"N_ENE",n_ene,max_ene) if (n_ene.gt.max_ene) then @@ -48,6 +51,13 @@ hamil_rep=index(controlcard,"HAMIL_REP").gt.0 write (iout,*) "Number of energy parameter sets",nparmset call multreadi(controlcard,"ISAMPL",isampl,nparmset,1) + do i=1,nparmset + if (isampl(i).eq.0) then + write (iout,*) "ERROR: isampl is 0 for parmset",i + call flush(iout) + stop + endif + enddo write (iout,*) "MaxSlice",MaxSlice call readi(controlcard,"NSLICE",nslice,1) call flush(iout) @@ -66,6 +76,8 @@ return1 endif indpdb=0 + energy_dec=(index(controlcard,'ENERGY_DEC').gt.0) + rmsrgymap = (index(controlcard,'RMSRGYMAP').gt.0) if (index(controlcard,"CLASSIFY").gt.0) indpdb=1 call reada(controlcard,"DELTA",delta,1.0d-2) call readi(controlcard,"EINICHECK",einicheck,2) @@ -73,7 +85,33 @@ call reada(controlcard,"DELTRGY",deltrgy,5.0d-2) call readi(controlcard,"RESCALE",rescale_mode,1) check_conf=index(controlcard,"NO_CHECK_CONF").eq.0 + call readi(controlcard,'TORMODE',tor_mode,0) +C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then + write(iout,*) "torsional and valence angle mode",tor_mode call reada(controlcard,'DISTCHAINMAX',distchainmax,50.0d0) + call reada(controlcard,'BOXX',boxxsize,100.0d0) + call reada(controlcard,'BOXY',boxysize,100.0d0) + call reada(controlcard,'BOXZ',boxzsize,100.0d0) +c Cutoff range for interactions + call reada(controlcard,"R_CUT",r_cut,15.0d0) + call reada(controlcard,"LAMBDA",rlamb,0.3d0) + call reada(controlcard,"LIPTHICK",lipthick,0.0d0) + call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0) + if (lipthick.gt.0.0d0) then + bordliptop=(boxzsize+lipthick)/2.0 + bordlipbot=bordliptop-lipthick +C endif + if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) + & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE" + buflipbot=bordlipbot+lipbufthick + bufliptop=bordliptop-lipbufthick + if ((lipbufthick*2.0d0).gt.lipthick) + &write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF" + endif + write(iout,*) "bordliptop=",bordliptop + write(iout,*) "bordlipbot=",bordlipbot + write(iout,*) "bufliptop=",bufliptop + write(iout,*) "buflipbot=",buflipbot call readi(controlcard,'SYM',symetr,1) write (iout,*) "DISTCHAINMAX",distchainmax write (iout,*) "delta",delta @@ -90,7 +128,34 @@ zscfile=index(controlcard,"ZSCFILE").gt.0 with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0 write (iout,*) "with_dihed_constr ",with_dihed_constr + with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0 + write (iout,*) "with_theta_constr ",with_theta_constr + call readi(controlcard,'SHIELD',shield_mode,0) + write(iout,*) "shield_mode",shield_mode +C endif + call readi(controlcard,'TORMODE',tor_mode,0) + write(iout,*) "torsional and valence angle mode",tor_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 +C do i=1,ntyp +C long_r_sidechain(i)=vbldsc0(1,i) +C short_r_sidechain(i)=sigma0(i) +C enddo + buff_shield=1.0d0 + endif + call readi(controlcard,'CONSTR_DIST',constr_dist,0) + dyn_ss=index(controlcard,"DYN_SS").gt.0 + adaptive = index(controlcard,"ADAPTIVE").gt.0 return end c------------------------------------------------------------------------------ @@ -400,7 +465,7 @@ c------------------------------------------------------------------------------- external ilen,iroof double precision rmsdev,energia(0:max_ene),efree,eini,temp double precision prop(maxQ) - integer ntot_all(maxslice,0:maxprocs-1) + integer ntot_all(maxslice,0:maxprocs-1), maxslice_buff integer iparm,ib,iib,ir,nprop,nthr,npars double precision etot,time integer ixdrf,iret @@ -531,7 +596,13 @@ c DA scratchfile. #ifdef MPI c Check if everyone has the same number of conformations - call MPI_Allgather(stot(1),maxslice,MPI_INTEGER, + +c call MPI_ALLgather(MPI_IN_PLACE,stot(1),MPI_DATATYPE_NULL, +c & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR) + + maxslice_buff=maxslice + + call MPI_Allgather(stot(1),maxslice_buff,MPI_INTEGER, & ntot_all(1,0),maxslice,MPI_INTEGER,MPI_Comm_World,IERROR) lerr=.false. do i=0,nprocs-1 @@ -772,3 +843,156 @@ c------------------------------------------------------------------------------- iroof = ii return end +c------------------------------------------------------------------------------- + subroutine read_dist_constr + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.CONTROL' + include 'COMMON.CHAIN' + include 'COMMON.IOUNITS' + include 'COMMON.SBRIDGE' + integer ifrag_(2,100),ipair_(2,100) + double precision wfrag_(100),wpair_(100) + character*500 controlcard + logical normalize + 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,.true.) + call readi(controlcard,"NFRAG",nfrag_,0) + call readi(controlcard,"NPAIR",npair_,0) + call readi(controlcard,"NDIST",ndist_,0) + call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) + call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0) + call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0) + call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0) + call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0) + normalize = index(controlcard,"NORMALIZE").gt.0 +c write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_ +c write (iout,*) "IFRAG" +c do i=1,nfrag_ +c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) +c enddo +c write (iout,*) "IPAIR" +c do i=1,npair_ +c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i) +c enddo + call flush(iout) + do i=1,nfrag_ + if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup + if (ifrag_(2,i).gt.nstart_sup+nsup-1) + & ifrag_(2,i)=nstart_sup+nsup-1 +c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i) + call flush(iout) + if (wfrag_(i).gt.0.0d0) then + do j=ifrag_(1,i),ifrag_(2,i)-1 + do k=j+1,ifrag_(2,i) +c write (iout,*) "j",j," k",k + ddjk=dist(j,k) + if (constr_dist.eq.1) then + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i) + else if (constr_dist.eq.2) then + if (ddjk.le.dist_cut) then + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i) + endif + else + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + dhpb(nhpb)=ddjk + forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2) + endif + write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + enddo + enddo + endif + enddo + do i=1,npair_ + if (wpair_(i).gt.0.0d0) then + ii = ipair_(1,i) + jj = ipair_(2,i) + if (ii.gt.jj) then + itemp=ii + ii=jj + jj=itemp + endif + do j=ifrag_(1,ii),ifrag_(2,ii) + do k=ifrag_(1,jj),ifrag_(2,jj) + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k + forcon(nhpb)=wpair_(i) + dhpb(nhpb)=dist(j,k) + write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) + enddo + enddo + endif + enddo + write (iout,*) "Distance restraints as read from input" + do i=1,ndist_ + 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) +c fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1) + if (forcon(nhpb+1).le.0.0d0.or.fordepth(nhpb+1).le.0.0d0)cycle + nhpb=nhpb+1 + write (iout,'(a,4i5,2f8.2,2f10.5)') "+dist.restr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(nhpb),dhpb(nhpb), + & dhpb1(nhpb),forcon(nhpb),fordepth(nhpb) + if (ibecarb(i).gt.0) then + ihpb(i)=ihpb(i)+nres + jhpb(i)=jhpb(i)+nres + endif + else +C print *,"in else" + read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i), + & ibecarb(i),forcon(nhpb+1) + if (forcon(nhpb+1).gt.0.0d0) then + nhpb=nhpb+1 + 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 + write (iout,'(a,4i5,f8.2,f10.1)') "+dist.restr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),ibecarb(i),dhpb(nhpb),forcon(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)) + enddo + if (constr_dist.eq.11 .and. normalize) then + fordepthmax=fordepth(1) + do i=2,nhpb + if (fordepth(i).gt.fordepthmax) fordepthmax=fordepth(i) + enddo + do i=1,nhpb + fordepth(i)=fordepth(i)/fordepthmax + enddo + write (iout,'(/a/4a5,2a8,2a10)') + & "Normalized Lorenzian-like distance restraints", + & " Nr"," res1"," res2"," beta"," d1"," d2"," k"," V" + do i=1,nhpb + write (iout,'(4i5,2f8.2,2f10.5)')i,ihpb(i),jhpb(i),ibecarb(i), + & dhpb(i),dhpb1(i),forcon(i),fordepth(i) + enddo + endif + write (iout,*) + call hpb_partition + call flush(iout) + return + end