X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Freadrtns_CSA.F;h=6984aeb5a837e74434529de48da582fae5b91acf;hb=b9025a366733ac0290debbd52ede3527d29c9ed5;hp=65f398069b3a28af5a9c889de709981d3f92296f;hpb=478a9d9a1c99eb3f4bc4ca676ff3162bdd01d633;p=unres.git diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 65f3980..6984aeb 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -139,7 +139,7 @@ C Set up the time limit (caution! The time must be input in minutes!) call readi(controlcard,'MAXGEN',maxgen,10000) call readi(controlcard,'MAXOVERLAP',maxoverlap,1000) call readi(controlcard,"KDIAG",kdiag,0) - call readi(controlcard,"RESCALE_MODE",rescale_mode,1) + call readi(controlcard,"RESCALE_MODE",rescale_mode,2) if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) & write (iout,*) "RESCALE_MODE",rescale_mode split_ene=index(controlcard,'SPLIT_ENE').gt.0 @@ -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) @@ -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, @@ -638,7 +639,7 @@ C 12/1/95 Added weight for the multi-body term WCORR & 'General scaling factor of SC-p interactions:',scalscp endif r0_corr=cutoff_corr-delt_corr - do i=1,20 + do i=1,ntyp aad(i,1)=scalscp*aad(i,1) aad(i,2)=scalscp*aad(i,2) bad(i,1)=scalscp*bad(i,1) @@ -679,13 +680,41 @@ 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,"DTRISS",dtriss,1.0D0) + call reada(weightcard,"ATRISS",atriss,0.3D0) + call reada(weightcard,"BTRISS",btriss,0.02D0) + call reada(weightcard,"CTRISS",ctriss,1.0D0) + 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 + ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain + 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 +726,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 @@ -717,7 +748,7 @@ c print *,'Finished reading pdb data' maxsi=1000 do i=2,nres-1 iti=itype(i) - if (iti.ne.10 .and. itype(i).ne.21) then + if (iti.ne.10 .and. itype(i).ne.ntyp1) then nsi=0 fail=.true. do while (fail.and.nsi.le.maxsi) @@ -749,8 +780,8 @@ C Assign initial virtual bond lengths vbld_inv(i)=vblinv enddo do i=2,nres-1 - vbld(i+nres)=dsc(itype(i)) - vbld_inv(i+nres)=dsc_inv(itype(i)) + vbld(i+nres)=dsc(iabs(itype(i))) + vbld_inv(i+nres)=dsc_inv(iabs(itype(i))) c write (iout,*) "i",i," itype",itype(i), c & " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres) enddo @@ -759,15 +790,15 @@ c print *,nres c print '(20i4)',(itype(i),i=1,nres) do i=1,nres #ifdef PROCOR - if (itype(i).eq.21 .or. itype(i+1).eq.21) then + if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then #else - if (itype(i).eq.21) then + if (itype(i).eq.ntyp1) then #endif itel(i)=0 #ifdef PROCOR - else if (itype(i+1).ne.20) then + else if (iabs(itype(i+1)).ne.20) then #else - else if (itype(i).ne.20) then + else if (iabs(itype(i)).ne.20) then #endif itel(i)=1 else @@ -824,8 +855,8 @@ C 8/13/98 Set limits to generating the dihedral angles #endif nct=nres cd print *,'NNT=',NNT,' NCT=',NCT - if (itype(1).eq.21) nnt=2 - if (itype(nres).eq.21) nct=nct-1 + if (itype(1).eq.ntyp1) nnt=2 + if (itype(nres).eq.ntyp1) nct=nct-1 if (pdbref) then if(me.eq.king.or..not.out1file) & write (iout,'(a,i3)') 'nsup=',nsup @@ -942,7 +973,7 @@ C initial geometry. enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then do j=1,3 dc(j,i+nres)=c(j,i+nres)-c(j,i) dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres) @@ -974,6 +1005,7 @@ C initial geometry. enddo do i=2,nres-1 omeg(i)=-120d0*deg2rad + if (itype(i).le.0) omeg(i)=-omeg(i) enddo else if(me.eq.king.or..not.out1file) @@ -1054,18 +1086,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) @@ -1122,17 +1171,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) @@ -1143,7 +1194,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 @@ -1211,7 +1263,7 @@ c enddo enddo do i=nnt,nct - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then 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) @@ -1291,7 +1343,7 @@ C Set up variable list. nvar=ntheta+nphi nside=0 do i=2,nres-1 - if (itype(i).ne.10 .and. itype(i).ne.21) then + if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then nside=nside+1 ialph(i,1)=nvar+nside ialph(nside,2)=i @@ -1650,6 +1702,9 @@ c---------------------------------------------------------------------------- call readi(minimcard,'MINFUN',minfun,maxmin) call reada(minimcard,'TOLF',tolf,1.0D-2) call reada(minimcard,'RTOLF',rtolf,1.0D-4) + print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1) + print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1) + print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1) write (iout,'(/80(1h*)/20x,a/80(1h*))') & 'Options in energy minimization:' write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') @@ -1912,12 +1967,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 @@ -2079,6 +2144,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) @@ -2220,11 +2287,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 @@ -2402,7 +2469,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.