X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fio_base.F90;h=f66a9b5d093b8cf48225cd7583e972e35917e300;hb=ae2f9918e1e5f88a78cc32b1d3a7ff90dede1207;hp=198fb1f0bd37b4b0f3e238047bf02a0cc4158262;hpb=0f5304c0d44e3ec4d6538de7d274c4b1a1930049;p=unres4.git

diff --git a/source/unres/io_base.F90 b/source/unres/io_base.F90
index 198fb1f..f66a9b5 100644
--- a/source/unres/io_base.F90
+++ b/source/unres/io_base.F90
@@ -435,7 +435,8 @@
       use control_data, only:out1file
       use MPI_data
       character*320 afmcard
-      integer i
+      real, dimension(3) ::cbeg,cend
+      integer i,j
       print *, "wchodze"
       call card_concat(afmcard,.true.)
       call readi(afmcard,"BEG",afmbeg,0)
@@ -443,15 +444,62 @@
       call reada(afmcard,"FORCE",forceAFMconst,0.0d0)
       call reada(afmcard,"VEL",velAFMconst,0.0d0)
       print *,'FORCE=' ,forceAFMconst
+      cbeg=0.0d0
+      cend=0.0d0
+      if (afmbeg.eq.-1) then
+        read(inp,*) nbegafmmat,(afmbegcentr(i),i=1,nbegafmmat)
+        do i=1,nbegafmmat
+         do j=1,3
+          cbeg(j)=cbeg(j)+c(j,afmbegcentr(i))/nbegafmmat
+         enddo
+        enddo
+      else
+      do j=1,3
+        cbeg(j)=c(j,afmend)
+      enddo
+      endif
+      if (afmend.eq.-1) then
+        read(inp,*) nendafmmat,(afmendcentr(i),i=1,nendafmmat)
+        do i=1,nendafmmat
+         do j=1,3
+          cend(j)=cend(j)+c(j,afmendcentr(i))/nendafmmat
+         enddo
+        enddo
+      else
+        cend(j)=c(j,afmend)
+      endif
 !------ NOW PROPERTIES FOR AFM
        distafminit=0.0d0
        do i=1,3
-        distafminit=(c(i,afmend)-c(i,afmbeg))**2+distafminit
+        distafminit=(cend(i)-cbeg(i))**2+distafminit
        enddo
         distafminit=dsqrt(distafminit)
         print *,'initdist',distafminit
       return
       end subroutine read_afminp
+!C------------------------------------------------------------
+      subroutine read_afmnano
+      use geometry_data
+      use energy_data
+      use control_data, only:out1file
+      use MPI_data
+      integer :: i
+      real*8 :: cm
+      read(inp,*,err=11) inanomove,(inanotab(i),i=1,inanomove),forcenanoconst 
+       cm=0.0d0
+       do i=1,inanomove
+         cm=cm+c(3,inanotab(i))
+       enddo
+       cm=cm/inanomove
+        distnanoinit=cm-tubecenter(3)
+      return
+11    write(iout,*)&
+      "error in afmnano",&
+       ", number of center, their index and forceconstance"
+      stop
+      return
+      end subroutine read_afmnano
+
 !-----------------------------------------------------------------------------
       subroutine read_dist_constr
       use MPI_data
@@ -475,9 +523,9 @@
       character(len=640) :: controlcard
 
 !el local variables
-      integer :: i,k,j,ddjk,ii,jj,itemp
+      integer :: i,k,j,ii,jj,itemp,mnumkk,mnumjj,k1,j1
       integer :: nfrag_,npair_,ndist_
-      real(kind=8) :: dist_cut
+      real(kind=8) :: dist_cut,ddjk
 
 !      write (iout,*) "Calling read_dist_constr"
 !      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
@@ -502,11 +550,11 @@
       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
-!      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
+      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
 !      write (iout,*) "IFRAG"
-!      do i=1,nfrag_
-!        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
-!      enddo
+      do i=1,nfrag_
+        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
+      enddo
 !      write (iout,*) "IPAIR"
 !      do i=1,npair_
 !        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
@@ -518,15 +566,25 @@
       
       call flush(iout)
       do i=1,nfrag_
-        if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
-        if (ifrag_(2,i).gt.nstart_sup+nsup-1) &
-          ifrag_(2,i)=nstart_sup+nsup-1
+!        if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
+!        if (ifrag_(2,i).gt.nstart_sup+nsup-1) &
+!          ifrag_(2,i)=nstart_sup+nsup-1
 !        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
         call flush(iout)
         if (wfrag_(i).gt.0.0d0) then
         do j=ifrag_(1,i),ifrag_(2,i)-1
           do k=j+1,ifrag_(2,i)
-!            write (iout,*) "j",j," k",k
+          write (iout,*) "j",j," k",k,nres
+          j1=j
+          k1=k
+          if (j.gt.nres) j1=j-nres
+          if (k.gt.nres) k1=k-nres
+          mnumkk=molnum(k1)
+          mnumjj=molnum(j1)
+          
+          if ((itype(k1,mnumkk).gt.ntyp_molec(mnumkk)).or.&
+              (itype(j1,mnumjj).gt.ntyp_molec(mnumjj))) cycle
+            write (iout,*) "j",j," k",k,itype(k1,mnumkk),itype(j1,mnumjj)
             ddjk=dist(j,k)
             if (constr_dist.eq.1) then
               nhpb=nhpb+1
@@ -1081,8 +1139,10 @@
       
       iatom=0
       ichain=1
-      ires=0
-      do i=nnt,nct
+      ires=0 
+      write(iout,*) "TUTUT"
+      do i=nnt,nct 
+        write(iout,*), "coord",c(1,i),c(2,i),c(3,i)
         iti=itype(i,molnum(i))
         print *,i,molnum(i)
         if (molnum(i+1).eq.0) then
@@ -1115,12 +1175,12 @@
         xi=c(1,i)
         yi=c(2,i)
         zi=c(3,i)
-        xi=dmod(xi,boxxsize)
-        if (xi.lt.0.0d0) xi=xi+boxxsize
-        yi=dmod(yi,boxysize)
-        if (yi.lt.0.0d0) yi=yi+boxysize
-        zi=dmod(zi,boxzsize)
-        if (zi.lt.0.0d0) zi=zi+boxzsize
+!        xi=dmod(xi,boxxsize)
+!        if (xi.lt.0.0d0) xi=xi+boxxsize
+!        yi=dmod(yi,boxysize)
+!        if (yi.lt.0.0d0) yi=yi+boxysize
+!        zi=dmod(zi,boxzsize)
+!        if (zi.lt.0.0d0) zi=zi+boxzsize
         write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
            ires,xi,yi,zi,vtot(i)
         else 
@@ -1552,7 +1612,7 @@
        write(iout,*) "Explicit LIPID not yet implemented"
        stop
       else if (molecule.eq.5) then
-      do i=1,ntyp1_molec(molecule)
+      do i=-1,ntyp1_molec(molecule)
         print *,restyp(i,molecule)
         print *,i,restyp(i,molecule)(1:2),nam(1:2)
         if (ucase(nam(1:2)).eq.restyp(i,molecule)(1:2)) then