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)
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
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
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)
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
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
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
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