X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Freadrtns_CSA.F;h=1cda5eac13362dfce355bd88dea59760e7c1eb10;hb=c65d296f3cdc480d3737eae335a510332b4f0dc8;hp=7ac066950029bae9fec1057eebf2d193201ebc7e;hpb=f5cea66c6a23f13718f902416e527a955693d213;p=unres.git diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 7ac0669..1cda5ea 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -295,11 +295,11 @@ cd endif do i=1,nrep iremd_m_total=iremd_m_total+remd_m(i) enddo - write (iout,*) 'Total number of replicas ',iremd_m_total + write (iout,*) 'Total number of replicas ',iremd_m_total + endif endif - endif if(me.eq.king.or..not.out1file) - & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup " + & write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup " return end c-------------------------------------------------------------------------- @@ -345,6 +345,7 @@ C rest = index(controlcard,"REST").gt.0 tbf = index(controlcard,"TBF").gt.0 usampl = index(controlcard,"USAMPL").gt.0 + mdpdb = index(controlcard,"MDPDB").gt.0 call reada(controlcard,"T_BATH",t_bath,300.0d0) call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1) @@ -538,7 +539,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 @@ -550,54 +551,54 @@ C C Body C C Read weights of the subsequent energy terms. - call card_concat(weightcard) - call reada(weightcard,'WLONG',wlong,1.0D0) - call reada(weightcard,'WSC',wsc,wlong) - call reada(weightcard,'WSCP',wscp,wlong) - call reada(weightcard,'WELEC',welec,1.0D0) - call reada(weightcard,'WVDWPP',wvdwpp,welec) - call reada(weightcard,'WEL_LOC',wel_loc,1.0D0) - call reada(weightcard,'WCORR4',wcorr4,0.0D0) - call reada(weightcard,'WCORR5',wcorr5,0.0D0) - call reada(weightcard,'WCORR6',wcorr6,0.0D0) - call reada(weightcard,'WTURN3',wturn3,1.0D0) - call reada(weightcard,'WTURN4',wturn4,1.0D0) - call reada(weightcard,'WTURN6',wturn6,1.0D0) - call reada(weightcard,'WSCCOR',wsccor,1.0D0) - call reada(weightcard,'WSTRAIN',wstrain,1.0D0) - call reada(weightcard,'WBOND',wbond,1.0D0) - call reada(weightcard,'WTOR',wtor,1.0D0) - 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,'SCAL14',scal14,0.4D0) - call reada(weightcard,'SCALSCP',scalscp,1.0d0) - call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) - call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) - call reada(weightcard,'TEMP0',temp0,300.0d0) - if (index(weightcard,'SOFT').gt.0) ipot=6 + call card_concat(weightcard) + call reada(weightcard,'WLONG',wlong,1.0D0) + call reada(weightcard,'WSC',wsc,wlong) + call reada(weightcard,'WSCP',wscp,wlong) + call reada(weightcard,'WELEC',welec,1.0D0) + call reada(weightcard,'WVDWPP',wvdwpp,welec) + call reada(weightcard,'WEL_LOC',wel_loc,1.0D0) + call reada(weightcard,'WCORR4',wcorr4,0.0D0) + call reada(weightcard,'WCORR5',wcorr5,0.0D0) + call reada(weightcard,'WCORR6',wcorr6,0.0D0) + call reada(weightcard,'WTURN3',wturn3,1.0D0) + call reada(weightcard,'WTURN4',wturn4,1.0D0) + call reada(weightcard,'WTURN6',wturn6,1.0D0) + call reada(weightcard,'WSCCOR',wsccor,1.0D0) + call reada(weightcard,'WSTRAIN',wstrain,1.0D0) + call reada(weightcard,'WBOND',wbond,1.0D0) + call reada(weightcard,'WTOR',wtor,1.0D0) + 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,'SCAL14',scal14,0.4D0) + call reada(weightcard,'SCALSCP',scalscp,1.0d0) + call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) + call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) + call reada(weightcard,'TEMP0',temp0,300.0d0) + if (index(weightcard,'SOFT').gt.0) ipot=6 C 12/1/95 Added weight for the multi-body term WCORR - call reada(weightcard,'WCORRH',wcorr,1.0D0) - if (wcorr4.gt.0.0d0) wcorr=wcorr4 - weights(1)=wsc - weights(2)=wscp - weights(3)=welec - weights(4)=wcorr - weights(5)=wcorr5 - weights(6)=wcorr6 - weights(7)=wel_loc - weights(8)=wturn3 - weights(9)=wturn4 - weights(10)=wturn6 - weights(11)=wang - weights(12)=wscloc - weights(13)=wtor - weights(14)=wtor_d - weights(15)=wstrain - weights(16)=wvdwpp - weights(17)=wbond - weights(18)=scal14 - weights(21)=wsccor + call reada(weightcard,'WCORRH',wcorr,1.0D0) + if (wcorr4.gt.0.0d0) wcorr=wcorr4 + weights(1)=wsc + weights(2)=wscp + weights(3)=welec + weights(4)=wcorr + weights(5)=wcorr5 + weights(6)=wcorr6 + weights(7)=wel_loc + weights(8)=wturn3 + weights(9)=wturn4 + weights(10)=wturn6 + weights(11)=wang + weights(12)=wscloc + weights(13)=wtor + weights(14)=wtor_d + weights(15)=wstrain + weights(16)=wvdwpp + weights(17)=wbond + weights(18)=scal14 + weights(21)=wsccor 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, @@ -679,13 +680,49 @@ 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. + enddo + do i=1,maxres-1 + do j=i+1,maxres + dyn_ssbond_ij(i,j)=1.0d300 + enddo + enddo + call reada(weightcard,"HT",Ht,0.0D0) + if (dyn_ss) then + ss_depth=ebr/wsc-0.25*eps(1,1) + Ht=Ht/wsc-0.25*eps(1,1) + akcm=akcm*wstrain/wsc + akth=akth*wstrain/wsc + akct=akct*wstrain/wsc + v1ss=v1ss*wstrain/wsc + v2ss=v2ss*wstrain/wsc + v3ss=v3ss*wstrain/wsc + else + 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 write (iout,*) "Parameters of the SS-bond potential:" write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth, & " AKCT",akct write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss - write (iout,*) "EBR",ebr -c print *,'indpdb=',indpdb,' pdbref=',pdbref + write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth + write (iout,*)" HT",Ht + print *,'indpdb=',indpdb,' pdbref=',pdbref endif if (indpdb.gt.0 .or. pdbref) then read(inp,'(a)') pdbfile @@ -697,9 +734,11 @@ c print *,'indpdb=',indpdb,' pdbref=',pdbref 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 @@ -789,21 +828,23 @@ 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) @@ -896,7 +937,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 @@ -916,7 +959,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 @@ -1055,18 +1098,35 @@ C Generate distance constraints, if the PDB structure is to be regularized. write (iout,'(/a,i3,a)') & 'The chain contains',ns,' disulfide-bridging cysteines.' write (iout,'(20i4)') (iss(i),i=1,ns) + if (dyn_ss) then + write(iout,*)"Running with dynamic disulfide-bond formation" + else write (iout,'(/a/)') 'Pre-formed links are:' do i=1,nss i1=ihpb(i)-nres 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) @@ -1082,6 +1142,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-------------------------------------------------------------------------- @@ -1123,17 +1184,19 @@ C Read information about disulfide bridges. include 'COMMON.SETUP' C Read bridging residues. read (inp,*) ns,(iss(i),i=1,ns) -c print *,'ns=',ns + print *,'ns=',ns if(me.eq.king.or..not.out1file) & write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns) 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) @@ -1144,7 +1207,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 @@ -1916,12 +1980,22 @@ C Get parameter filenames and open the parameter files. open (itordp,file=tordname,status='old',readonly) call getenv_loc('SCCORPAR',sccorname) open (isccor,file=sccorname,status='old',readonly) +#ifndef CRYST_THETA + call getenv_loc('THETPARPDB',thetname_pdb) + print *,"thetname_pdb ",thetname_pdb + open (ithep_pdb,file=thetname_pdb,status='old',action='read') + print *,ithep_pdb," opened" +#endif call getenv_loc('FOURIER',fouriername) open (ifourier,file=fouriername,status='old',readonly) call getenv_loc('ELEPAR',elename) open (ielep,file=elename,status='old',readonly) call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',readonly) +#ifndef CRYST_SC + call getenv_loc('ROTPARPDB',rotname_pdb) + open (irotam_pdb,file=rotname_pdb,status='old',action='read') +#endif #endif #ifndef OLDSCP C @@ -2083,6 +2157,8 @@ C Write file names & thetname(:ilen(thetname)) write (iout,*) "Rotamer parameter file : ", & rotname(:ilen(rotname)) + write (iout,*) "Thetpdb parameter file : ", + & thetname_pdb(:ilen(thetname_pdb)) write (iout,*) "Threading database : ", & patname(:ilen(patname)) if (lentmp.ne.0) @@ -2190,7 +2266,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) @@ -2224,11 +2301,11 @@ c write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(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 + nhpb=nhpb+1 + ihpb(nhpb)=j + jhpb(nhpb)=k dhpb(nhpb)=ddjk - forcon(nhpb)=wfrag_(i) + forcon(nhpb)=wfrag_(i) else if (constr_dist.eq.2) then if (ddjk.le.dist_cut) then nhpb=nhpb+1 @@ -2284,11 +2361,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 ", @@ -2297,7 +2393,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 @@ -2406,7 +2502,7 @@ C #endif if (OKRandom) then c r1 = prng_next(me) - r1=ran_number(0.0D0,1.0D0) + r1=ran_number(0.0D0,1.0D0) if(me.eq.king) & write (iout,*) 'ran_num',r1 if (r1.lt.0.0d0) OKRandom=.false.