if (refstr) then
call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,4f12.3,i5,$)')&
- itime,totT,EK,potE,totE,&
+ itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
potEcomp(23),me
format1="a133"
else
!C print *,'A CHUJ',potEcomp(23)
write (line1,'(i10,f15.2,7f12.3,i5,$)') &
- itime,totT,EK,potE,totE,&
+ itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
kinetic_T,t_bath,gyrate(),&
potEcomp(23),me
format1="a114"
call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2, &
f9.3,i5,$)') &
- itime,totT,EK,potE,totE,&
+ itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
distance,potEcomp(23),me
format1="a133"
else
!C print *,'A CHUJ',potEcomp(23)
write (line1,'(i10,f15.2,8f12.3,i5,$)')&
- itime,totT,EK,potE,totE, &
+ itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51), &
kinetic_T,t_bath,gyrate(),&
distance,potEcomp(23),me
format1="a114"
endif
+ else if (velnanoconst.ne.0 ) then
+ if (refstr) then
+ call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
+ write (line1,'(i10,f15.2,3f12.3,f7.2,2f6.3,f12.3,f10.1,2f8.2, &
+ f9.3,i5,$)') &
+ itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
+ rms,frac,frac_nn,kinetic_T,t_bath,gyrate(),&
+ vecsim,vectrue,me
+ format1="a133"
+!C print *,"CHUJOWO"
+ else
+!C print *,'A CHUJ',potEcomp(23)
+ write (line1,'(i10,f15.2,8f12.3,i5,$)')&
+ itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51), &
+ kinetic_T,t_bath,gyrate(),&
+ vecsim,vectrue,me
+ format1="a114"
+ endif
else
if (refstr) then
call rms_nac_nnc(rms,frac,frac_nn,co,.false.)
write (line1,'(i10,f15.2,3f12.3,f7.2,4f6.3,3f12.3,i5,$)') &
- itime,totT,EK,potE,totE,&
+ itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
rms,frac,frac_nn,co,amax,kinetic_T,t_bath,gyrate(),me
format1="a133"
else
write (line1,'(i10,f15.2,7f12.3,i5,$)') &
- itime,totT,EK,potE,totE,&
+ itime,totT,EK,potE+potEcomp(51),totE+potEcomp(51),&
amax,kinetic_T,t_bath,gyrate(),me
format1="a114"
endif
use MPI_data
use compare, only:seq_comp,contact
use control
+ implicit none
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
#ifdef MPI
! integer ilen
!el external ilen
!el local varables
- integer :: i,j,l,k,kkk,ii,i1,i2,it1,it2
+ integer :: i,j,l,k,kkk,ii,i1,i2,it1,it2,mnum
real(kind=8),dimension(3,maxres2+2) :: c_alloc
real(kind=8),dimension(3,0:maxres2) :: dc_alloc
call reada(weightcard,'WCATCAT',wcatcat,0.0d0)
call reada(weightcard,'WCATPROT',wcatprot,0.0d0)
call reada(weightcard,'WCATNUCL',wcatnucl,0.0d0)
+ call reada(weightcard,'WCATTRAN',wcat_tran,0.0d0)
call reada(weightcard,'WPEPBASE',wpepbase,1.0d0)
call reada(weightcard,'WSCPHO',wscpho,0.0d0)
call reada(weightcard,'WPEPPHO',wpeppho,0.0d0)
weights(48)=wscpho
weights(49)=wpeppho
weights(50)=wcatnucl
-
+ weights(56)=wcat_tran
if(me.eq.king.or..not.out1file) &
write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,&
wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,&
! print *,'Finished reading pdb data'
if(me.eq.king.or..not.out1file) &
- write (iout,'(a,i3,a,i3)')'nsup=',nsup,&
+ write (iout,'(a,i7,a,i7)')'nsup=',nsup,&
' nstart_sup=',nstart_sup !,"ergwergewrgae"
!el if(.not.allocated(itype_pdb))
allocate(itype_pdb(nres))
endif
if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
do i=2,nres
+ mnum=molnum(i)
if (molnum(i).eq.1) then
vbld(i)=vbl
vbld_inv(i)=vblinv
else
- vbld(i)=7.0
- vbld_inv(i)=1.0/7.0
+ write(iout,*) "typy",itype(i,mnum),ntyp1_molec(mnum),i
+ vbld(i)=6.0
+ vbld_inv(i)=1.0/6.0
+ if ((itype(i,mnum).eq.ntyp1_molec(mnum)).or.&
+ (itype(i-1,mnum).eq.ntyp1_molec(mnum))) then
+ vbld(i)=1.9
+ vbld_inv(i)=1.0/1.9
+ endif
endif
enddo
do i=2,nres-1
+ mnum=molnum(i)
if (molnum(i).eq.1) then
! print *, "molnum",molnum(i),itype(i,molnum(i)),nres,i
vbld(i+nres)=dsc(iabs(itype(i,molnum(i))))
vbld_inv(i+nres)=dsc_inv(iabs(itype(i,molnum(i))))
else
+ write(iout,*) "typy2",itype(i,mnum),ntyp1_molec(mnum),i
+ if (itype(i,mnum).eq.ntyp1_molec(mnum)) then
+ vbld_inv(i+nres)=1.0
+ vbld(i+nres)=0.0
+ else
vbld(i+nres)=vbldsc0_nucl(1,iabs(itype(i,molnum(i))))
vbld_inv(i+nres)=1.0/vbldsc0_nucl(1,iabs(itype(i,molnum(i))))
+ endif
endif
! write (iout,*) "i",i," itype",itype(i,1),
! & " dsc",dsc(itype(i,1))," vbld",vbld(i),vbld(i+nres)
! print *,"tu dochodze"
! znamy nres oraz nss można zaalokowac potrzebne tablice
call alloc_geo_arrays
+ print *,"after alloc_geo_arrays"
call alloc_ener_arrays
+ print *,"after alloc_ener_arrays"
!--------------------------------
! 8/13/98 Set limits to generating the dihedral angles
do i=1,nres
endif
#endif
nct=nres
+ allocate(ireschain(nres))
+ ireschain=0
+ write(iout,*),"before seq2chains",ireschain
+ call seq2chains
+ write(iout,*) "after seq2chains",nchain
+ allocate ( chain_border1(2,nchain))
+ chain_border1(1,1)=1
+ chain_border1(2,1)=chain_border(2,1)+1
+ do i=2,nchain-1
+ chain_border1(1,i)=chain_border(1,i)-1
+ chain_border1(2,i)=chain_border(2,i)+1
+ enddo
+ if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
+ chain_border1(2,nchain)=nres
+ write(iout,*) "nres",nres," nchain",nchain
+ do i=1,nchain
+ write(iout,*)"chain",i,chain_length(i),chain_border(1,i),&
+ chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
+ enddo
+ allocate(tabpermchain(50,5040))
+ call chain_symmetry(npermchain,tabpermchain)
print *,'NNT=',NNT,' NCT=',NCT
if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
if (pdbref) then
if(me.eq.king.or..not.out1file) &
- write (iout,'(a,i3)') 'nsup=',nsup
+ write (iout,'(a,i7)') 'nsup=',nsup
nstart_seq=nnt
if (nsup.le.(nct-nnt+1)) then
do i=0,nct-nnt+1-nsup
'Error - sequences to be superposed do not match.'
endif
111 continue
- if (nsup.eq.0) nsup=nct-nnt
+ if (nsup.eq.0) nsup=nct-nnt+1
if (nstart_sup.eq.0) nstart_sup=nnt
if (nstart_seq.eq.0) nstart_seq=nnt
if(me.eq.king.or..not.out1file) &
endif
if(me.eq.king.or..not.out1file)then
write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
- write (iout,*) 'IZ_SC=',iz_sc
+ write (iout,*) 'IZ_SC=',iz_sc, 'NCT=',nct
endif
!----------------------
call init_int_table
if (constr_dist.gt.0) call read_dist_constr
write (iout,*) "After read_dist_constr nhpb",nhpb
if ((AFMlog.gt.0).or.(selfguide.gt.0)) call read_afminp
+ if (velnanoconst.ne.0) call read_afmnano
call hpb_partition
if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4 &
write (iout,'(a)') 'Extended chain initial geometry.'
do i=3,nres
theta(i)=90d0*deg2rad
+ if (molnum(i).eq.2) theta(i)=160d0*deg2rad
enddo
do i=4,nres
phi(i)=180d0*deg2rad
+ if (molnum(i).eq.2) phi(i)=30d0*deg2rad
enddo
do i=2,nres-1
alph(i)=110d0*deg2rad
+ if (molnum(i).eq.2) alph(i)=30d0*deg2rad
enddo
do i=2,nres-1
omeg(i)=-120d0*deg2rad
+ if (molnum(i).eq.2) omeg(i)=60d0*deg2rad
if (itype(i,1).le.0) omeg(i)=-omeg(i)
enddo
call chainbuild
! enddo
ii_in_use=0
if(.not.allocated(pdbfiles_chomo)) allocate(pdbfiles_chomo(constr_homology))
- if(.not.allocated(chomo)) allocate(chomo(3,nres,constr_homology))
+ if(.not.allocated(chomo)) allocate(chomo(3,2*nres+2,constr_homology))
do k=1,constr_homology
end subroutine read_constr_homology
!-----------------------------------------------------------------------------
subroutine read_klapaucjusz
-! implicit none
+ use energy_data
+ implicit none
! include 'DIMENSIONS'
!#ifdef MPI
! include 'mpif.h'
return
10 stop "Error in fragment file"
end subroutine read_klapaucjusz
+
+!-----------------------------------------------------------
+ subroutine seq2chains
+!c
+!c Split the total UNRES sequence, which has dummy residues separating
+!c the chains, into separate chains. The length of chain ichain is
+!c contained in chain_length(ichain), the first and last non-dummy
+!c residues are in chain_border(1,ichain) and chain_border(2,ichain),
+!c respectively. The lengths pertain to non-dummy residues only.
+!c
+! implicit none
+! include 'DIMENSIONS'
+ use energy_data, only:molnum,nchain,chain_length,ireschain
+ implicit none
+! integer ireschain(nres)
+ integer ii,ichain,i,j,mnum
+ logical new_chain
+ print *,"in seq2"
+ ichain=1
+ new_chain=.true.
+ if (.not.allocated(chain_length)) allocate(chain_length(50))
+ if (.not.allocated(chain_border)) allocate(chain_border(2,50))
+
+ chain_length(ichain)=0
+ ii=1
+ do while (ii.lt.nres)
+ write(iout,*) "in seq2chains",ii,new_chain
+ mnum=molnum(ii)
+ if (itype(ii,mnum).eq.ntyp1_molec(mnum)) then
+ if (.not.new_chain) then
+ new_chain=.true.
+ chain_border(2,ichain)=ii-1
+ ichain=ichain+1
+ chain_border(1,ichain)=ii+1
+ chain_length(ichain)=0
+ endif
+ else
+ if (new_chain) then
+ chain_border(1,ichain)=ii
+ new_chain=.false.
+ endif
+ chain_length(ichain)=chain_length(ichain)+1
+ endif
+ ii=ii+1
+ enddo
+ if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) then
+ ii=ii-1
+ else
+ chain_length(ichain)=chain_length(ichain)+1
+ endif
+ if (chain_length(ichain).gt.0) then
+ chain_border(2,ichain)=ii
+ nchain=ichain
+ else
+ nchain=ichain-1
+ endif
+ ireschain=0
+ do i=1,nchain
+ do j=chain_border(1,i),chain_border(2,i)
+ ireschain(j)=i
+ enddo
+ enddo
+ return
+ end subroutine
+!---------------------------------------------------------------------
+ subroutine chain_symmetry(npermchain,tabpermchain)
+!c
+!c Determine chain symmetry. nperm is the number of permutations and
+!c tabperchain contains the allowed permutations of the chains.
+!c
+! implicit none
+! include "DIMENSIONS"
+! include "COMMON.IOUNITS"
+ implicit none
+ integer itemp(50),&
+ npermchain,tabpermchain(50,5040),&
+ tabperm(50,5040),mapchain(50),&
+ iequiv(50,nres),iflag(nres)
+ integer i,j,k,l,ii,nchain_group,nequiv(50),iieq,&
+ nperm,npermc,ind,mnum
+ if (nchain.eq.1) then
+ npermchain=1
+ tabpermchain(1,1)=1
+!c print*,"npermchain",npermchain," tabpermchain",tabpermchain(1,1)
+ return
+ endif
+!c
+!c Look for equivalent chains
+#ifdef DEBUG
+ write(iout,*) "nchain",nchain
+ do i=1,nchain
+ write(iout,*) "chain",i," from",chain_border(1,i),&
+ " to",chain_border(2,i)
+ write(iout,*)&
+ "sequence ",(itype(j,molnum(j)),j=chain_border(1,i),chain_border(2,i))
+ enddo
+#endif
+ do i=1,nchain
+ iflag(i)=0
+ enddo
+ nchain_group=0
+ do i=1,nchain
+ if (iflag(i).gt.0) cycle
+ iflag(i)=1
+ nchain_group=nchain_group+1
+ iieq=1
+ iequiv(iieq,nchain_group)=i
+ do j=i+1,nchain
+ if (iflag(j).gt.0.or.chain_length(i).ne.chain_length(j)) cycle
+!c k=0
+!c do while(k.lt.chain_length(i) .and.
+!c & itype(chain_border(1,i)+k).eq.itype(chain_border(1,j)+k))
+ do k=0,chain_length(i)-1
+!c k=k+1
+ mnum=molnum(k)
+ if (itype(chain_border(1,i)+k,mnum).ne.&
+ itype(chain_border(1,j)+k,mnum)) exit
+ enddo
+ if (k.lt.chain_length(i)) cycle
+ iflag(j)=1
+ iieq=iieq+1
+ iequiv(iieq,nchain_group)=j
+ enddo
+ nequiv(nchain_group)=iieq
+ enddo
+ write(iout,*) "Number of equivalent chain groups:",nchain_group
+ write(iout,*) "Equivalent chain groups"
+ do i=1,nchain_group
+ write(iout,*) "group",i," #members",nequiv(i)," chains",&
+ (iequiv(j,i),j=1,nequiv(i))
+ enddo
+ ind=0
+ do i=1,nchain_group
+ do j=1,nequiv(i)
+ ind=ind+1
+ mapchain(ind)=iequiv(j,i)
+ enddo
+ enddo
+ write (iout,*) "mapchain"
+ do i=1,nchain
+ write (iout,*) i,mapchain(i)
+ enddo
+ ii=0
+ do i=1,nchain_group
+ call permut(nequiv(i),nperm,tabperm)
+ if (ii.eq.0) then
+ ii=nequiv(i)
+ npermchain=nperm
+ do j=1,nperm
+ do k=1,ii
+ tabpermchain(k,j)=iequiv(tabperm(k,j),i)
+ enddo
+ enddo
+ else
+ npermc=npermchain
+ npermchain=npermchain*nperm
+ ind=0
+ do k=1,nperm
+ do j=1,npermc
+ ind=ind+1
+ do l=1,ii
+ tabpermchain(l,ind)=tabpermchain(l,j)
+ enddo
+ do l=1,nequiv(i)
+ tabpermchain(ii+l,ind)=iequiv(tabperm(l,k),i)
+ enddo
+ enddo
+ enddo
+ ii=ii+nequiv(i)
+ endif
+ enddo
+ do i=1,npermchain
+ do j=1,nchain
+ itemp(mapchain(j))=tabpermchain(j,i)
+ enddo
+ do j=1,nchain
+ tabpermchain(j,i)=itemp(j)
+ enddo
+ enddo
+ write(iout,*) "Number of chain permutations",npermchain
+ write(iout,*) "Permutations"
+ do i=1,npermchain
+ write(iout,'(20i4)') (tabpermchain(j,i),j=1,nchain)
+ enddo
+ return
+ end
+!c---------------------------------------------------------------------
+ integer function tperm(i,iperm,tabpermchain)
+! implicit none
+! include 'DIMENSIONS'
+ integer i,iperm
+ integer tabpermchain(50,5040)
+ if (i.eq.0) then
+ tperm=0
+ else
+ tperm=tabpermchain(i,iperm)
+ endif
+ return
+ end function
+
!-----------------------------------------------------------------------------
end module io
projects / unres4.git / blobdiff