X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Freadrtns.F;h=10cde78274ba4d8d9e3ef5dc006b77eb1e84a9ec;hb=e123f68beb7954771018d0311468c68250b5f2d1;hp=8f70874b6f849099d7a6671df9fe624201167327;hpb=fbc0f12bd8022b2c766eb8c350e7bb0451816c1d;p=unres.git diff --git a/source/unres/src_MD/readrtns.F b/source/unres/src_MD/readrtns.F index 8f70874..10cde78 100644 --- a/source/unres/src_MD/readrtns.F +++ b/source/unres/src_MD/readrtns.F @@ -138,7 +138,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 @@ -860,12 +860,36 @@ 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) + 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 + 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 @@ -893,6 +917,9 @@ c print *,'Finished reading pdb data' call contact(.false.,ncont_ref,icont_ref,co) if (sideadd) then +C Following 2 lines for diagnostics; comment out if not needed + write (iout,*) "Before sideadd" + call intout if(me.eq.king.or..not.out1file) & write(iout,*)'Adding sidechains' maxsi=1000 @@ -909,7 +936,12 @@ c print *,'Finished reading pdb data' & i,' after ',nsi,' trials' endif enddo +C 10/03/12 Adam: Recalculate coordinates with new side chain positions + call chainbuild endif +C Following 2 lines for diagnostics; comment out if not needed + write (iout,*) "After sideadd" + call intout endif if (indpdb.eq.0) then C Read sequence if not taken from the pdb file. @@ -1095,10 +1127,10 @@ czscore call geom_to_var(nvar,coord_exp_zs(1,1)) c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup if (constr_dist.gt.0) then call read_dist_constr - call hpb_partition endif + if (nhpb.gt.0) call hpb_partition c write (iout,*) "After read_dist_constr nhpb",nhpb - call flush(iout) +c call flush(iout) 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 @@ -1217,6 +1249,21 @@ 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" + 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 + else write (iout,'(/a/)') 'Pre-formed links are:' do i=1,nss i1=ihpb(i)-nres @@ -1229,6 +1276,7 @@ C Generate distance constraints, if the PDB structure is to be regularized. & ebr,forcon(i) enddo write (iout,'(a)') + endif endif if (i2ndstr.gt.0) call secstrp2dihc c call geom_to_var(nvar,x) @@ -1292,10 +1340,12 @@ 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) @@ -1813,6 +1863,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)') @@ -2083,38 +2136,38 @@ C Get parameter filenames and open the parameter files. open (isidep,file=sidename,status='old') #else open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old', - & readonly) + &action='read') open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown') C open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown') C Get parameter filenames and open the parameter files. call getenv_loc('BONDPAR',bondname) - open (ibond,file=bondname,status='old',readonly) + open (ibond,file=bondname,status='old',action='read') call getenv_loc('THETPAR',thetname) - open (ithep,file=thetname,status='old',readonly) + open (ithep,file=thetname,status='old',action='read') #ifndef CRYST_THETA call getenv_loc('THETPARPDB',thetname_pdb) print *,"thetname_pdb ",thetname_pdb - open (ithep_pdb,file=thetname_pdb,status='old',readonly) + open (ithep_pdb,file=thetname_pdb,status='old',action='read') print *,ithep_pdb," opened" #endif call getenv_loc('ROTPAR',rotname) - open (irotam,file=rotname,status='old',readonly) + open (irotam,file=rotname,status='old',action='read') #ifndef CRYST_SC call getenv_loc('ROTPARPDB',rotname_pdb) - open (irotam_pdb,file=rotname_pdb,status='old',readonly) + open (irotam_pdb,file=rotname_pdb,status='old',action='read') #endif call getenv_loc('TORPAR',torname) - open (itorp,file=torname,status='old',readonly) + open (itorp,file=torname,status='old',action='read') call getenv_loc('TORDPAR',tordname) - open (itordp,file=tordname,status='old',readonly) + open (itordp,file=tordname,status='old',action='read') call getenv_loc('SCCORPAR',sccorname) - open (isccor,file=sccorname,status='old',readonly) + open (isccor,file=sccorname,status='old',action='read') call getenv_loc('FOURIER',fouriername) - open (ifourier,file=fouriername,status='old',readonly) + open (ifourier,file=fouriername,status='old',action='read') call getenv_loc('ELEPAR',elename) - open (ielep,file=elename,status='old',readonly) + open (ielep,file=elename,status='old',action='read') call getenv_loc('SIDEPAR',sidename) - open (isidep,file=sidename,status='old',readonly) + open (isidep,file=sidename,status='old',action='read') #endif #ifndef OLDSCP C @@ -2129,7 +2182,7 @@ C #elif (defined G77) open (iscpp,file=scpname,status='old') #else - open (iscpp,file=scpname,status='old',readonly) + open (iscpp,file=scpname,status='old',action='read') #endif #endif call getenv_loc('PATTERN',patname) @@ -2140,7 +2193,7 @@ C #elif (defined G77) open (icbase,file=patname,status='old') #else - open (icbase,file=patname,status='old',readonly) + open (icbase,file=patname,status='old',action='read') #endif #ifdef MPI C Open output file only for CG processes @@ -2195,6 +2248,7 @@ c print *,"Processor",myrank," fg_rank",fg_rank #if defined(AIX) || defined(PGI) if (me.eq.king .or. .not. out1file) & open(iout,file=outname,status='unknown') +c#define DEBUG #ifdef DEBUG if (fg_rank.gt.0) then write (liczba,'(i3.3)') myrank/nfgtasks @@ -2203,6 +2257,7 @@ c print *,"Processor",myrank," fg_rank",fg_rank & status='unknown') endif #endif +c#undef DEBUG if(me.eq.king) then open(igeom,file=intname,status='unknown',position='append') open(ipdb,file=pdbname,status='unknown') @@ -2214,14 +2269,17 @@ c1out open(iout,file=outname,status='unknown') #else if (me.eq.king .or. .not.out1file) & open(iout,file=outname,status='unknown') +c#define DEBUG #ifdef DEBUG if (fg_rank.gt.0) then + print "Processor",fg_rank," opening output file" write (liczba,'(i3.3)') myrank/nfgtasks write (ll,'(bz,i3.3)') fg_rank open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll, & status='unknown') endif #endif +c#undef DEBUG if(me.eq.king) then open(igeom,file=intname,status='unknown',access='append') open(ipdb,file=pdbname,status='unknown')