X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Freadrtns_CSA.F;h=8172868618480b7d3f201ed33bca18cc5b73c014;hb=e12c8db731e266c6e38b2f5883f3ae7c15bbfb98;hp=79840a4579fa4e63dd95ecc5e767a1d143bed3ff;hpb=aedd2f72960d9cec69fcc7f0b6fd6555c7e22312;p=unres.git diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 79840a4..8172868 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 @@ -179,6 +183,7 @@ C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then 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) @@ -262,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 @@ -409,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 @@ -531,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" @@ -600,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) @@ -641,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) @@ -750,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. @@ -767,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 @@ -780,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, @@ -1038,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 @@ -1111,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.' @@ -1504,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)) @@ -1514,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) @@ -1911,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 @@ -2127,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 @@ -2327,10 +2427,11 @@ c------------------------------------------------------------------------------- open(irest2,file=rest2name,status='unknown') read(irest2,*) totT,EK,potE,totE,t_bath totTafm=totT - do i=1,2*nres + do i=0,2*nres read(irest2,'(3e15.5)') (d_t(j,i),j=1,3) enddo - do i=1,2*nres +C WARING IN OUTER FIELD CRITICAL CORRECTION TO CHANGE i=0 + do i=0,2*nres read(irest2,'(3e15.5)') (dc(j,i),j=1,3) enddo if(usampl) then @@ -2431,6 +2532,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) @@ -2556,6 +2661,7 @@ C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) #endif enddo + 1712 continue call flush(iout) return end