introduction nanosphere part 1
[unres.git] / source / unres / src_MD-M / readrtns_CSA.F
index e40717f..ce4c59f 100644 (file)
@@ -147,6 +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"
+      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 
@@ -268,6 +270,7 @@ C      endif
       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)
@@ -422,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
@@ -544,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"
@@ -766,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.
@@ -783,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
@@ -1521,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))
@@ -1531,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)
@@ -1928,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
@@ -2450,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)
@@ -2575,6 +2660,7 @@ C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
 #endif
 
       enddo
+ 1712 continue
       call flush(iout)
       return
       end