implicit none
!-----------------------------------------------------------------------------
! Max. number of AA residues
- integer,parameter :: maxres=6500!1200
+ integer,parameter :: maxres=101000!1200
! Appr. max. number of interaction sites
integer,parameter :: maxres2=2*maxres
! parameter (maxres6=6*maxres)
!-----------------------------------------------------------------------------
subroutine read_x(kanal,*)
- use geometry_data
- use energy_data
- use geometry, only:int_from_cart1
+ use geometry_data
+ use energy_data
+ use geometry, only:int_from_cart1
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'COMMON.GEO'
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
!-----------------------------------------------------------------------------
subroutine pdbout(etot,tytul,iunit)
- use geometry_data, only: c,nres
+ use geometry_data, only: c,nres,boxxsize,boxysize,boxzsize
use energy_data
! use control
use compare_data
! include 'COMMON.MD'
!el character(len=50) :: tytul
character*(*) :: tytul
- character(len=1),dimension(10) :: chainid= (/'A','B','C','D','E','F','G','H','I','J'/)
+ character(len=1),dimension(58) :: chainid= (/'A','B','C','D','E','F','G','H','I','J',&
+ 'K','L','M','O','Q','P','R','S','T','U','W','X','Y','Z','a','b','c','d','e','f',&
+ 'g','h','i','j','k','l','m','n','o','p','r','s','t','u','w','x','y','z',&
+ '0','1','2','3','4','5','6','7','8','9'/)
#ifdef XDRFPDB
integer,dimension(maxres) :: ica !(maxres)
#else
integer :: iti1
!el local variables
integer :: j,iti,itj,itk,itl,i,iatom,ichain,ires,iunit
- real(kind=8) :: etot
+ real(kind=8) :: etot,xi,yi,zi
integer :: nres2
nres2=2*nres
#ifdef XDRFPDB
if(.not.allocated(molnum)) allocate(molnum(maxres))
- molnum(:)=1
+! molnum(:)=mnumi(:)
if(.not.allocated(vtot)) allocate(vtot(maxres*2)) !(maxres2)
+ if(.not.allocated(istype)) allocate(istype(maxres))
#else
if(.not.allocated(vtot)) allocate(vtot(nres2)) !(maxres2)
#endif
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
endif
if (iti.eq.ntyp1_molec(molnum(i))) then
ichain=ichain+1
+ if (ichain.gt.58) ichain=1
! ires=0
write (iunit,'(a)') 'TER'
else
if (istype(i).eq.0) istype(i)=1
write (iunit,40) iatom,sugartyp(istype(i)),restyp(iti,2), &
chainid(ichain),ires,(c(j,i),j=1,3),vtot(i)
+ elseif (molnum(i).eq.4) 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
+ write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
+ ires,xi,yi,zi,vtot(i)
+ elseif (molnum(i).eq.5) 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
+ write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
+ ires,xi,yi,zi,vtot(i)
else
+
write (iunit,60) iatom,restyp(iti,molnum(i)),chainid(ichain),&
ires,(c(j,i),j=1,3),vtot(i)
endif
enddo
write (iunit,'(a)') 'TER'
do i=nnt,nct-1
+ if (molnum(i).eq.5) cycle
if (itype(i,1).eq.ntyp1) cycle
if ((itype(i,1).eq.10 .and. itype(i+1,1).ne.ntyp1).or.(molnum(i).eq.5)) then
write (iunit,30) ica(i),ica(i+1)
do i=1,nres
iti=itype(i,1)
! print *,vbld(i),"vbld(i)",i
- write (iout,'(a3,i4,6f10.3)') restyp(iti,1),i,vbld(i),&
+ write (iout,'(a3,i6,6f10.3)') restyp(iti,1),i,vbld(i),&
rad2deg*theta(i),rad2deg*phi(i),vbld(nres+i),rad2deg*alph(i),&
rad2deg*omeg(i)
enddo
!-----------------------------------------------------------------------------
! permut.F
!-----------------------------------------------------------------------------
- subroutine permut(isym)
-
- use geometry_data, only: tabperm
-! implicit real*8 (a-h,o-z)
-! include 'DIMENSIONS'
-! include 'COMMON.LOCAL'
-! include 'COMMON.VAR'
-! include 'COMMON.CHAIN'
-! include 'COMMON.INTERACT'
-! include 'COMMON.IOUNITS'
-! include 'COMMON.GEO'
-! include 'COMMON.NAMES'
-! include 'COMMON.CONTROL'
-
- integer :: n,isym
+ subroutine permut(isym,nperm,tabperm)
+!c integer maxperm,maxsym
+!c parameter (maxperm=3628800)
+!c parameter (maxsym=10)
+! include "DIMENSIONS"
+ integer n,a,tabperm,nperm,kkk,i,isym
! logical nextp
-!el external nextp
- integer,dimension(isym) :: a
-! parameter(n=symetr)
-!el local variables
- integer :: kkk,i
-
+! external nextp
+ dimension a(isym),tabperm(50,5040)
n=isym
+ nperm=1
if (n.eq.1) then
tabperm(1,1)=1
return
endif
+ do i=2,n
+ nperm=nperm*i
+ enddo
kkk=0
do i=1,n
a(i)=i
enddo
- 10 print *,(a(i),i=1,n)
+ 10 continue
+!c print '(i3,2x,100i3)',kkk+1,(a(i),i=1,n)
kkk=kkk+1
do i=1,n
- tabperm(kkk,i)=a(i)
-! write (iout,*) "tututu", kkk
+ tabperm(i,kkk)=a(i)
enddo
if(nextp(n,a)) go to 10
return
- end subroutine permut
+ end subroutine
+
+
!-----------------------------------------------------------------------------
logical function nextp(n,a)
!-----------------------------------------------------------------------------
integer function rescode(iseq,nam,itype,molecule)
-! use io_base, only: ucase
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'COMMON.NAMES'
else if (molecule.eq.2) then
do i=1,ntyp1_molec(molecule)
print *,nam(1:1),restyp(i,molecule)(1:1)
+ print *,nam(2:2),"read",i
if (nam(2:2).eq.restyp(i,molecule)(1:1)) then
rescode=i
return
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
endif
return
end function sugarcode
+ integer function rescode_lip(res,ires)
+ character(len=3):: res
+ integer ires
+ rescode_lip=0
+ if (res.eq.'NC3') rescode_lip =4
+ if (res.eq.'NH3') rescode_lip =2
+ if (res.eq.'PO4') rescode_lip =3
+ if (res.eq.'GL1') rescode_lip =12
+ if (res.eq.'GL2') rescode_lip =12
+ if (res.eq.'C1A') rescode_lip =18
+ if (res.eq.'C2A') rescode_lip =18
+ if (res.eq.'C3A') rescode_lip =18
+ if (res.eq.'C4A') rescode_lip =18
+ if (res.eq.'C1B') rescode_lip =18
+ if (res.eq.'C2B') rescode_lip =18
+ if (res.eq.'C3B') rescode_lip =18
+ if (res.eq.'C4B') rescode_lip =18
+ if (res.eq.'D2A') rescode_lip =16
+ if (res.eq.'D2B') rescode_lip =16
+
+ if (rescode_lip.eq.0) write(iout,*) "UNKNOWN RESIDUE",ires,res
+ return
+ end function
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
end module io_base
projects / unres4.git / blobdiff