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=5ba2a6014fbd859202ed28885b66520ac09bb77d;hp=89ddeaca181c1bb74aa8bbc8bf291f8661211f37;hpb=c1827efa66e69c93c6e2b2e7420b06b430c3550a;p=unres.git diff --git a/source/unres/src_MD-M/readrtns_CSA.F b/source/unres/src_MD-M/readrtns_CSA.F index 89ddeac..ce4c59f 100644 --- a/source/unres/src_MD-M/readrtns_CSA.F +++ b/source/unres/src_MD-M/readrtns_CSA.F @@ -147,11 +147,8 @@ C constrains on theta angles WITH_THETA_CONSTR is the keyword print *,'AFMlog',AFMlog,selfguide,"KUPA" call readi(controlcard,'TUBEMOD',tubelog,0) write (iout,*) TUBElog,"TUBEMODE" - if (TUBElog.gt.0) then - call reada(controlcard,"XTUBE",tubecenter(1),0.0d0) - call reada(controlcard,"YTUBE",tubecenter(2),0.0d0) - call reada(controlcard,"RTUBE",tubeR0,0.0d0) - endif + 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 @@ -270,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 @@ -417,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 @@ -539,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" @@ -761,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. @@ -778,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 @@ -1516,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)) @@ -1526,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) @@ -1923,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 @@ -2445,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) @@ -2570,6 +2660,7 @@ C dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb)) #endif enddo + 1712 continue call flush(iout) return end