X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Freadrtns_CSA.F;h=ce4c59f25307b98a454799b1417f9357c5b7b82b;hb=d2ee52eab91e6055d9cd36a42447313d73b00ea1;hp=bdc7c9d029ba2373cdf7906f92e7a49d380e0ddb;hpb=2095e643aafa2fb192e066e69e5ca3b1a59c2727;p=unres.git diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index bdc7c9d..ce4c59f 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -104,7 +104,7 @@ 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 + call reada(controlcard,'TIMLIM',timlim,2800.0D0) ! default 16 hours unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes call reada(controlcard,'RMSDBC',rmsdbc,3.0D0) @@ -145,6 +145,10 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword AFMlog=(index(controlcard,'AFM')) selfguide=(index(controlcard,'SELFGUIDE')) print *,'AFMlog',AFMlog,selfguide,"KUPA" + call readi(controlcard,'TUBEMOD',tubelog,0) + write (iout,*) TUBElog,"TUBEMODE" + call readi(controlcard,'GENCONSTR',genconstr,0) +C write (iout,*) TUBElog,"TUBEMODE" 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 @@ -155,6 +159,9 @@ C 2 reseved for further possible developement 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,'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 readi(controlcard,'MAXGEN',maxgen,10000) call readi(controlcard,'MAXOVERLAP',maxoverlap,1000) call readi(controlcard,"KDIAG",kdiag,0) @@ -176,6 +183,7 @@ C endif else icheckgrad=3 endif + call reada(controlcard,'DELTA',aincr,1.0d-7) elseif (index(controlcard,'THREAD').gt.0) then modecalc=2 call readi(controlcard,'THREAD',nthread,0) @@ -247,7 +255,7 @@ c Cutoff range for interactions bordliptop=(boxzsize+lipthick)/2.0 bordlipbot=bordliptop-lipthick C endif - if ((bordliptop.gt.boxzsize).or.(borlipbot.lt.0.0)) + if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) & write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE" buflipbot=bordlipbot+lipbufthick bufliptop=bordliptop-lipbufthick @@ -259,6 +267,17 @@ C endif write(iout,*) "bufliptop=",bufliptop write(iout,*) "buflipbot=",buflipbot write (iout,*) "SHIELD MODE",shield_mode + if (TUBElog.gt.0) then + call reada(controlcard,"XTUBE",tubecenter(1),0.0d0) + call reada(controlcard,"YTUBE",tubecenter(2),0.0d0) + call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0) + call reada(controlcard,"RTUBE",tubeR0,0.0d0) + call reada(controlcard,"TUBETOP",bordtubetop,boxzsize) + call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0) + call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0) + buftubebot=bordtubebot+tubebufthick + buftubetop=bordtubetop-tubebufthick + endif if (shield_mode.gt.0) then pi=3.141592d0 C VSolvSphere the volume of solving sphere @@ -406,6 +425,14 @@ c call reada(controlcard,"R_CUT",r_cut,2.0d0) c call reada(controlcard,"LAMBDA",rlamb,0.3d0) rest = index(controlcard,"REST").gt.0 tbf = index(controlcard,"TBF").gt.0 + tnp = index(controlcard,"NOSEPOINCARE99").gt.0 + tnp1 = index(controlcard,"NOSEPOINCARE01").gt.0 + tnh = index(controlcard,"NOSEHOOVER96").gt.0 + if (RESPA.and.tnh)then + xiresp = index(controlcard,"XIRESP").gt.0 + endif + write(iout,*) "tnh",tnh + call reada(controlcard,"Q_NP",Q_np,0.1d0) usampl = index(controlcard,"USAMPL").gt.0 mdpdb = index(controlcard,"MDPDB").gt.0 @@ -528,6 +555,43 @@ c Calculate friction coefficients and bounds of stochastic forces & "Velocities will be reset at random every",count_reset_vel, & " steps" endif + else if (tnp .or. tnp1 .or. tnh) then + if (tnp .or. tnp1) then + write (iout,'(a)') "Nose-Poincare bath calculation" + if (tnp) write (iout,'(a)') + & "J.Comput.Phys. 151 114 (1999) S.D.Bond B.J.Leimkuhler B.B.Laird" + if (tnp1) write (iout,'(a)') "JPSJ 70 75 (2001) S. Nose" + else + write (iout,'(a)') "Nose-Hoover bath calculation" + write (iout,'(a)') "Mol.Phys. 87 1117 (1996) Martyna et al." + nresn=1 + nyosh=1 + nnos=1 + do i=1,nnos + qmass(i)=Q_np + xlogs(i)=1.0 + vlogs(i)=0.0 + enddo + do i=1,nyosh + WDTI(i) = 1.0*d_time/nresn + WDTI2(i)=WDTI(i)/2 + WDTI4(i)=WDTI(i)/4 + WDTI8(i)=WDTI(i)/8 + enddo + if (RESPA) then + if(xiresp) then + write (iout,'(a)') "NVT-XI-RESPA algorithm" + else + write (iout,'(a)') "NVT-XO-RESPA algorithm" + endif + do i=1,nyosh + WDTIi(i) = 1.0*d_time/nresn/ntime_split + WDTIi2(i)=WDTIi(i)/2 + WDTIi4(i)=WDTIi(i)/4 + WDTIi8(i)=WDTIi(i)/8 + enddo + endif + endif else if(me.eq.king.or..not.out1file) & write (iout,'(a31)') "Microcanonical mode calculation" @@ -597,6 +661,7 @@ C include 'COMMON.BOUNDS' include 'COMMON.MD' include 'COMMON.SETUP' + include 'COMMON.SHIELD' character*4 sequence(maxres) integer rescode double precision x(maxvar) @@ -638,7 +703,9 @@ C Read weights of the subsequent energy terms. call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0) call reada(weightcard,'DELT_CORR',delt_corr,0.5d0) call reada(weightcard,'TEMP0',temp0,300.0d0) + call reada(weightcard,'WSHIELD',wshield,1.0d0) call reada(weightcard,'WLT',wliptran,0.0D0) + call reada(weightcard,'WTUBE',wtube,1.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) @@ -747,10 +814,16 @@ C 12/1/95 Added weight for the multi-body term WCORR call reada(weightcard,"BTRISS",btriss,0.021D0) call reada(weightcard,"CTRISS",ctriss,1.001D0) call reada(weightcard,"DTRISS",dtriss,1.001D0) + call reada(weightcard,"LIPSCALE",lipscale,1.3D0) write (iout,*) "ATRISS=", atriss write (iout,*) "BTRISS=", btriss write (iout,*) "CTRISS=", ctriss write (iout,*) "DTRISS=", dtriss + if (shield_mode.gt.0) then + lipscale=0.0d0 + write (iout,*) "liscale not used in shield mode" + endif + write (iout,*) "lipid scaling factor", lipscale dyn_ss=(index(weightcard,'DYN_SS').gt.0) do i=1,maxres dyn_ss_mask(i)=.false. @@ -764,12 +837,12 @@ C 12/1/95 Added weight for the multi-body term WCORR 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 + akcm=akcm/wsc + akth=akth/wsc + akct=akct/wsc + v1ss=v1ss/wsc + v2ss=v2ss/wsc + v3ss=v3ss/wsc else if (wstrain.ne.0.0) then ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain @@ -777,7 +850,7 @@ C 12/1/95 Added weight for the multi-body term WCORR ss_depth=0.0 endif endif - + write (iout,*) "wshield,", wshield if(me.eq.king.or..not.out1file) then write (iout,*) "Parameters of the SS-bond potential:" write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth, @@ -1035,7 +1108,7 @@ c---------------------- call MPI_Finalize(MPI_COMM_WORLD,IERROR) stop 'Error reading reference structure' #endif - 39 call chainbuild + 39 call chainbuild_extconf call setup_var czscore call geom_to_var(nvar,coord_exp_zs(1,1)) nstart_sup=nnt @@ -1108,6 +1181,7 @@ C initial geometry. return else call read_angles(inp,*36) + call chainbuild_extconf endif goto 37 36 write (iout,'(a)') 'Error reading angle file.' @@ -1132,6 +1206,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.' @@ -1500,6 +1575,31 @@ cd print *,'gen_dist_constr: nnt=',nnt,' nct=',nct cd write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct, cd & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq, cd & ' nsup',nsup + if (constr_dist.eq.11) then + do i=nstart_sup,nstart_sup+nsup-1 + do j=i+2,nstart_sup+nsup-1 + distance=dist(i,j) + if (distance.le.15.0) then + nhpb=nhpb+1 + ihpb(nhpb)=i+nstart_seq-nstart_sup + jhpb(nhpb)=j+nstart_seq-nstart_sup + forcon(nhpb)=sqrt(0.04*distance) + fordepth(nhpb)=sqrt(40.0/distance) + dhpb(nhpb)=distance-0.1d0 + dhpb1(nhpb)=distance+0.1d0 + +#ifdef MPI + if (.not.out1file .or. me.eq.king) + & write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) +#else + write (iout,'(a,3i5,f8.2,f10.2)') "+dist.constr ", + & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb) +#endif + endif + enddo + enddo + else do i=nstart_sup,nstart_sup+nsup-1 cd write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)), cd & ' seq_pdb', restyp(itype_pdb(i)) @@ -1510,7 +1610,8 @@ cd & ' seq_pdb', restyp(itype_pdb(i)) forcon(nhpb)=weidis dhpb(nhpb)=dist(i,j) enddo - enddo + enddo + endif cd write (iout,'(a)') 'Distance constraints:' cd do i=nss+1,nhpb cd ii=ihpb(i) @@ -1907,6 +2008,7 @@ c---------------------------------------------------------------------------- iread=index(rekord,lancuch(:ilen(lancuch))//"=") if (iread.eq.0) return iread=iread+ilen(lancuch)+1 + read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim) 10 return end @@ -1999,12 +2101,18 @@ C Get parameter filenames and open the parameter files. open (itorp,file=torname,status='old',readonly,shared) call getenv_loc('TORDPAR',tordname) open (itordp,file=tordname,status='old',readonly,shared) + call getenv_loc('TORKCC',torkccname) + open (itorkcc,file=torkccname,status='old',readonly,shared) + call getenv_loc('THETKCC',thetkccname) + open (ithetkcc,file=thetkccname,status='old',readonly,shared) call getenv_loc('FOURIER',fouriername) open (ifourier,file=fouriername,status='old',readonly,shared) call getenv_loc('ELEPAR',elename) open (ielep,file=elename,status='old',readonly,shared) call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',readonly,shared) + call getenv_loc('LIPTRANPAR',liptranname) + open (iliptranpar,file=liptranname,status='old',readonly,shared) #elif (defined CRAY) || (defined AIX) open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old', & action='read') @@ -2027,6 +2135,10 @@ c print *,"Processor",myrank," opened file IROTAM" c print *,"Processor",myrank," opened file ITORP" call getenv_loc('TORDPAR',tordname) open (itordp,file=tordname,status='old',action='read') + call getenv_loc('TORKCC',torkccname) + open (itorkcc,file=torkccname,status='old',action='read') + call getenv_loc('THETKCC',thetkccname) + open (ithetkcc,file=thetkccname,status='old',action='read') c print *,"Processor",myrank," opened file ITORDP" call getenv_loc('SCCORPAR',sccorname) open (isccor,file=sccorname,status='old',action='read') @@ -2039,6 +2151,8 @@ c print *,"Processor",myrank," opened file IFOURIER" c print *,"Processor",myrank," opened file IELEP" call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old',action='read') + call getenv_loc('LIPTRANPAR',liptranname) + open (iliptranpar,file=liptranname,status='old',action='read') c print *,"Processor",myrank," opened file ISIDEP" c print *,"Processor",myrank," opened parameter files" #elif (defined G77) @@ -2056,6 +2170,10 @@ C Get parameter filenames and open the parameter files. open (itorp,file=torname,status='old') call getenv_loc('TORDPAR',tordname) open (itordp,file=tordname,status='old') + call getenv_loc('TORKCC',torkccname) + open (itorkcc,file=torkccname,status='old') + call getenv_loc('THETKCC',thetkccname) + open (ithetkcc,file=thetkccname,status='old') call getenv_loc('SCCORPAR',sccorname) open (isccor,file=sccorname,status='old') call getenv_loc('FOURIER',fouriername) @@ -2064,6 +2182,8 @@ C Get parameter filenames and open the parameter files. open (ielep,file=elename,status='old') call getenv_loc('SIDEPAR',sidename) open (isidep,file=sidename,status='old') + call getenv_loc('LIPTRANPAR',liptranname) + open (iliptranpar,file=liptranname,status='old') #else open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old', & readonly) @@ -2080,6 +2200,10 @@ C Get parameter filenames and open the parameter files. open (itorp,file=torname,status='old',readonly) call getenv_loc('TORDPAR',tordname) open (itordp,file=tordname,status='old',readonly) + call getenv_loc('TORKCC',torkccname) + open (itorkcc,file=torkccname,status='old',readonly) + call getenv_loc('THETKCC',thetkccname) + open (ithetkcc,file=thetkccname,status='old',readonly) call getenv_loc('SCCORPAR',sccorname) open (isccor,file=sccorname,status='old',readonly) #ifndef CRYST_THETA @@ -2101,6 +2225,8 @@ C Get parameter filenames and open the parameter files. open (irotam_pdb,file=rotname_pdb,status='old',action='read') #endif #endif + call getenv_loc('TUBEPAR',tubename) + open (itube,file=tubename,status='old',readonly) #ifndef OLDSCP C C 8/9/01 In the newest version SCp interaction constants are read from a file @@ -2405,6 +2531,10 @@ c------------------------------------------------------------------------------- write (iout,*) "Calling read_dist_constr" c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup c call flush(iout) + if ((genconstr.gt.0).and.(constr_dist.eq.11)) then + call gen_dist_constr + go to 1712 + endif call card_concat(controlcard) call readi(controlcard,"NFRAG",nfrag_,0) call readi(controlcard,"NPAIR",npair_,0) @@ -2530,6 +2660,7 @@ C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) #endif enddo + 1712 continue call flush(iout) return end