use io_base
use geometry_data
use geometry
+ use control_data, only:maxterm_sccor
implicit none
!-----------------------------------------------------------------------------
! Max. number of residue types and parameters in expressions for
! parameter (maxtor=4,maxterm=10)
!-----------------------------------------------------------------------------
! Max number of torsional terms in SCCOR
- integer,parameter :: maxterm_sccor=6
+!el integer,parameter :: maxterm_sccor=6
!-----------------------------------------------------------------------------
character(len=1),dimension(:),allocatable :: secstruc !(maxres)
!-----------------------------------------------------------------------------
!
!-----------------------------------------------------------------------------
contains
-#ifndef WHAM_RUN
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
!-----------------------------------------------------------------------------
! bank.F io_csa
!-----------------------------------------------------------------------------
endif
enddo
ndih_constr=ii
- deallocate(secstruc)
+! deallocate(secstruc)
return
100 continue
write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
- deallocate(secstruc)
+! deallocate(secstruc)
return
110 continue
write(iout,'(A20)')'Error reading FTORS'
- deallocate(secstruc)
+! deallocate(secstruc)
return
end subroutine secstrp2dihc
!-----------------------------------------------------------------------------
' in position',i4)
return
end subroutine read_secstr_pred
-#endif
+!#endif
!-----------------------------------------------------------------------------
! parmread.F
!-----------------------------------------------------------------------------
use geometry_data
use energy_data
- use control_data, only:maxtor,maxterm
+ use control_data, only:maxterm !,maxtor
use MD_data
use MPI_data
!el use map_data
real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
- res1
+ res1,epsijlip
integer :: ichir1,ichir2
! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
!el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
!el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
-
!
! For printing parameters after they are read set the following in the UNRES
! C-shell script:
allocate(isc(ntyp+1)) !(ntyp+1)
allocate(restok(ntyp+1)) !(ntyp+1)
allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
+ allocate(long_r_sidechain(ntyp))
+ allocate(short_r_sidechain(ntyp))
+ dsc(:)=0.0d0
+ dsc_inv(:)=0.0d0
#ifdef CRYST_BOND
read (ibond,*) vbldp0,akp,mp,ip,pstok
endif
enddo
#else
- read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
+ read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp,ip,pstok
do i=1,ntyp
read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
j=1,nbondterm(i)),msc(i),isc(i),restok(i)
allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
- do i=-ntyp,ntyp
- a0thet(i)=0.0D0
- do j=1,2
- do ichir1=-1,1
- do ichir2=-1,1
- athet(j,i,ichir1,ichir2)=0.0D0
- bthet(j,i,ichir1,ichir2)=0.0D0
- enddo
- enddo
- enddo
- do j=0,3
- polthet(j,i)=0.0D0
- enddo
- do j=1,3
- gthet(j,i)=0.0D0
- enddo
- theta0(i)=0.0D0
- sig0(i)=0.0D0
- sigc0(i)=0.0D0
- enddo
+
+ a0thet(:)=0.0D0
+ athet(:,:,:,:)=0.0D0
+ bthet(:,:,:,:)=0.0D0
+ polthet(:,:)=0.0D0
+ gthet(:,:)=0.0D0
+ theta0(:)=0.0D0
+ sig0(:)=0.0D0
+ sigc0(:)=0.0D0
+ allocate(liptranene(ntyp))
+!C reading lipid parameters
+ write (iout,*) "iliptranpar",iliptranpar
+ call flush(iout)
+ read(iliptranpar,*) pepliptran
+ print *,pepliptran
+ do i=1,ntyp
+ read(iliptranpar,*) liptranene(i)
+ print *,liptranene(i)
+ enddo
+ close(iliptranpar)
#ifdef CRYST_THETA
!
!----------------------------------------------------
allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
- allocate(aa0thet(-maxthetyp1:maxthetyp1,&
- -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+ allocate(aa0thet(-nthetyp-1:nthetyp+1,&
+ -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
!(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
- allocate(aathet(ntheterm,-maxthetyp1:maxthetyp1,&
- -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+ allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
+ -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
!(maxtheterm,-maxthetyp1:maxthetyp1,&
! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
- allocate(bbthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
- -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
- allocate(ccthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
- -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
- allocate(ddthet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
- -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
- allocate(eethet(nsingle,ntheterm2,-maxthetyp1:maxthetyp1,&
- -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+ allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
+ -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
+ allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
+ -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
+ allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
+ -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
+ allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
+ -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
!(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
- allocate(ffthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,&
- -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
- allocate(ggthet(ndouble,ndouble,ntheterm3,-maxthetyp1:maxthetyp1,&
- -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+ allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
+ -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
+ allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
+ -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
!(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
do i=-ntyp1,-1
ithetyp(i)=-ithetyp(-i)
enddo
- do iblock=1,2
- do i=-maxthetyp,maxthetyp
- do j=-maxthetyp,maxthetyp
- do k=-maxthetyp,maxthetyp
- aa0thet(i,j,k,iblock)=0.0d0
- do l=1,ntheterm
- aathet(l,i,j,k,iblock)=0.0d0
- enddo
- do l=1,ntheterm2
- do m=1,nsingle
- bbthet(m,l,i,j,k,iblock)=0.0d0
- ccthet(m,l,i,j,k,iblock)=0.0d0
- ddthet(m,l,i,j,k,iblock)=0.0d0
- eethet(m,l,i,j,k,iblock)=0.0d0
- enddo
- enddo
- do l=1,ntheterm3
- do m=1,ndouble
- do mm=1,ndouble
- ffthet(mm,m,l,i,j,k,iblock)=0.0d0
- ggthet(mm,m,l,i,j,k,iblock)=0.0d0
- enddo
- enddo
- enddo
- enddo
- enddo
- enddo
- enddo
+
+ aa0thet(:,:,:,:)=0.0d0
+ aathet(:,:,:,:,:)=0.0d0
+ bbthet(:,:,:,:,:,:)=0.0d0
+ ccthet(:,:,:,:,:,:)=0.0d0
+ ddthet(:,:,:,:,:,:)=0.0d0
+ eethet(:,:,:,:,:,:)=0.0d0
+ ffthet(:,:,:,:,:,:,:)=0.0d0
+ ggthet(:,:,:,:,:,:,:)=0.0d0
+
! VAR:iblock means terminally blocking group 1=non-proline 2=proline
do iblock=1,2
! VAR:ntethtyp is type of theta potentials type currently 0=glycine
allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
- do i=1,ntyp
- do j=1,maxlob
- bsc(j,i)=0.0D0
- nlob(i)=0
- enddo
- enddo
- nlob(ntyp1)=0
- dsc(ntyp1)=0.0D0
-
- do i=-ntyp,ntyp
- do j=1,maxlob
- do k=1,3
- censc(k,j,i)=0.0D0
- enddo
- do k=1,3
- do l=1,3
- gaussc(l,k,j,i)=0.0D0
- enddo
- enddo
- enddo
- enddo
-
+ bsc(:,:)=0.0D0
+ nlob(:)=0
+ nlob(:)=0
+ dsc(:)=0.0D0
+ censc(:,:,:)=0.0D0
+ gaussc(:,:,:,:)=0.0D0
+
#ifdef CRYST_SC
!
! Read the parameters of the probability distribution/energy expression
allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
!el---------------------------
- do iblock=1,2
- do i=-ntortyp,ntortyp
- do j=-ntortyp,ntortyp
- nterm(i,j,iblock)=0
- nlor(i,j,iblock)=0
- enddo
- enddo
- enddo
+ nterm(:,:,:)=0
+ nlor(:,:,:)=0
!el---------------------------
read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
do i=-ntyp,-1
itortyp(i)=-itortyp(-i)
enddo
- write (iout,*) 'ntortyp',ntortyp
+ itortyp(ntyp1)=ntortyp
+ itortyp(-ntyp1)=-ntortyp
do iblock=1,2
+ write (iout,*) 'ntortyp',ntortyp
do i=0,ntortyp-1
do j=-ntortyp+1,ntortyp-1
read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
.or. t3.ne.toronelet(k)) then
- write (iout,*) "Error in double torsional parameter file",&
+ write (iout,*) "Error in double torsional parameter file",&
i,j,k,t1,t2,t3
#ifdef MPI
call MPI_Finalize(Ierror)
write (iout,*) "Coefficients of the cumulants"
endif
read (ifourier,*) nloctyp
-write(iout,*) "nloctyp",nloctyp
+!write(iout,*) "nloctyp",nloctyp
!el from module energy-------
allocate(b1(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
allocate(b2(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
allocate(ee(2,2,-nloctyp-1:nloctyp+1))
allocate(ctilde(2,2,-nloctyp-1:nloctyp+1))
allocate(dtilde(2,2,-nloctyp-1:nloctyp+1)) !(2,2,-maxtor:maxtor)
+! el
+ b1(1,:)=0.0d0
+ b1(2,:)=0.0d0
!--------------------------------
do i=0,nloctyp-1
read (ifourier,*,end=115,err=115)
read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
if (lprint) then
- write (iout,*) 'Type',i
- write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
+ write (iout,*) 'Type',i
+ write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
endif
B1(1,i) = b(3)
B1(2,i) = b(5)
enddo
enddo
endif
-
!
! Read electrostatic-interaction parameters
!
allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
- do i=1,ntyp
- do j=1,ntyp
- augm(i,j)=0.0D0
- enddo
- chip(i)=0.0D0
- alp(i)=0.0D0
- sigma0(i)=0.0D0
- sigii(i)=0.0D0
- rr0(i)=0.0D0
- enddo
+ allocate(epslip(ntyp,ntyp))
+ augm(:,:)=0.0D0
+ chip(:)=0.0D0
+ alp(:)=0.0D0
+ sigma0(:)=0.0D0
+ sigii(:)=0.0D0
+ rr0(:)=0.0D0
+
!--------------------------------
read (isidep,*,end=117,err=117) ipot,expon
if(me.eq.king) &
write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
', exponents are ',expon,2*expon
- goto (10,20,30,30,40) ipot
+! goto (10,20,30,30,40) ipot
+ select case(ipot)
!----------------------- LJ potential ---------------------------------
- 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
- (sigma0(i),i=1,ntyp)
- if (lprint) then
- write (iout,'(/a/)') 'Parameters of the LJ potential:'
- write (iout,'(a/)') 'The epsilon array:'
- call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
- write (iout,'(/a)') 'One-body parameters:'
- write (iout,'(a,4x,a)') 'residue','sigma'
- write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
- endif
- goto 50
+ case (1)
+! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+ read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+ (sigma0(i),i=1,ntyp)
+ if (lprint) then
+ write (iout,'(/a/)') 'Parameters of the LJ potential:'
+ write (iout,'(a/)') 'The epsilon array:'
+ call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
+ write (iout,'(/a)') 'One-body parameters:'
+ write (iout,'(a,4x,a)') 'residue','sigma'
+ write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
+ endif
+! goto 50
!----------------------- LJK potential --------------------------------
- 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
- (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
- if (lprint) then
- write (iout,'(/a/)') 'Parameters of the LJK potential:'
- write (iout,'(a/)') 'The epsilon array:'
- call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
- write (iout,'(/a)') 'One-body parameters:'
- write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
- write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
- i=1,ntyp)
- endif
- goto 50
+ case(2)
+! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+ read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+ (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
+ if (lprint) then
+ write (iout,'(/a/)') 'Parameters of the LJK potential:'
+ write (iout,'(a/)') 'The epsilon array:'
+ call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
+ write (iout,'(/a)') 'One-body parameters:'
+ write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
+ write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),&
+ i=1,ntyp)
+ endif
+! goto 50
!---------------------- GB or BP potential -----------------------------
- 30 do i=1,ntyp
- read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
- enddo
- read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
- read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
- read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
- read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
-! For the GB potential convert sigma'**2 into chi'
- if (ipot.eq.4) then
- do i=1,ntyp
- chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
+ case(3:4)
+! 30 do i=1,ntyp
+ do i=1,ntyp
+ read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
enddo
- endif
- if (lprint) then
- write (iout,'(/a/)') 'Parameters of the BP potential:'
- write (iout,'(a/)') 'The epsilon array:'
- call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
- write (iout,'(/a)') 'One-body parameters:'
- write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
- ' chip ',' alph '
- write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
- chip(i),alp(i),i=1,ntyp)
- endif
- goto 50
+ read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
+ read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
+ read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
+ read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
+ do i=1,ntyp
+ read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
+ enddo
+
+! For the GB potential convert sigma'**2 into chi'
+ if (ipot.eq.4) then
+ do i=1,ntyp
+ chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
+ enddo
+ endif
+ if (lprint) then
+ write (iout,'(/a/)') 'Parameters of the BP potential:'
+ write (iout,'(a/)') 'The epsilon array:'
+ call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
+ write (iout,'(/a)') 'One-body parameters:'
+ write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
+ ' chip ',' alph '
+ write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),&
+ chip(i),alp(i),i=1,ntyp)
+ endif
+! goto 50
!--------------------- GBV potential -----------------------------------
- 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
- (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
- (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
- if (lprint) then
- write (iout,'(/a/)') 'Parameters of the GBV potential:'
- write (iout,'(a/)') 'The epsilon array:'
- call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
- write (iout,'(/a)') 'One-body parameters:'
- write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
- 's||/s_|_^2',' chip ',' alph '
- write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
- sigii(i),chip(i),alp(i),i=1,ntyp)
- endif
- 50 continue
+ case(5)
+! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+ read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
+ (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
+ (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
+ if (lprint) then
+ write (iout,'(/a/)') 'Parameters of the GBV potential:'
+ write (iout,'(a/)') 'The epsilon array:'
+ call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
+ write (iout,'(/a)') 'One-body parameters:'
+ write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
+ 's||/s_|_^2',' chip ',' alph '
+ write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),&
+ sigii(i),chip(i),alp(i),i=1,ntyp)
+ endif
+ case default
+ write(iout,*)"Wrong ipot"
+! 50 continue
+ end select
+ continue
close (isidep)
!-----------------------------------------------------------------------
! Calculate the "working" parameters of SC interactions.
!el from module energy - COMMON.INTERACT-------
- allocate(aa(ntyp1,ntyp1),bb(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
+ allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
+ allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
allocate(sigma(0:ntyp1,0:ntyp1),r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
- do i=1,ntyp1
- do j=1,ntyp1
- aa(i,j)=0.0D0
- bb(i,j)=0.0D0
- chi(i,j)=0.0D0
- sigma(i,j)=0.0D0
- r0(i,j)=0.0D0
- enddo
- enddo
+ aa_aq(:,:)=0.0D0
+ bb_aq(:,:)=0.0D0
+ aa_lip(:,:)=0.0D0
+ bb_lip(:,:)=0.0D0
+ chi(:,:)=0.0D0
+ sigma(:,:)=0.0D0
+ r0(:,:)=0.0D0
+
!--------------------------------
do i=2,ntyp
do j=1,i-1
- eps(i,j)=eps(j,i)
+ eps(i,j)=eps(j,i)
+ epslip(i,j)=epslip(j,i)
enddo
enddo
do i=1,ntyp
epsij=eps(i,j)
sigeps=dsign(1.0D0,epsij)
epsij=dabs(epsij)
- aa(i,j)=epsij*rrij*rrij
- bb(i,j)=-sigeps*epsij*rrij
- aa(j,i)=aa(i,j)
- bb(j,i)=bb(i,j)
+ aa_aq(i,j)=epsij*rrij*rrij
+ bb_aq(i,j)=-sigeps*epsij*rrij
+ aa_aq(j,i)=aa_aq(i,j)
+ bb_aq(j,i)=bb_aq(i,j)
+ epsijlip=epslip(i,j)
+ sigeps=dsign(1.0D0,epsijlip)
+ epsijlip=dabs(epsijlip)
+ aa_lip(i,j)=epsijlip*rrij*rrij
+ bb_lip(i,j)=-sigeps*epsijlip*rrij
+ aa_lip(j,i)=aa_lip(i,j)
+ bb_lip(j,i)=bb_lip(i,j)
+!C write(iout,*) aa_lip
if (ipot.gt.2) then
sigt1sq=sigma0(i)**2
sigt2sq=sigma0(j)**2
endif
if (lprint) then
write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
- restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),&
+ restyp(i),restyp(j),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
endif
enddo
enddo
allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
- do i=1,ntyp
- do j=1,2
- bad(i,j)=0.0D0
- enddo
- enddo
+ bad(:,:)=0.0D0
#ifdef OLDSCP
!
stop
return
end subroutine parmread
-!-----------------------------------------------------------------------------
-! permut.F
-!-----------------------------------------------------------------------------
- subroutine permut(isym)
-
- use geometry_data, only: tabperm
-! use energy_data
-! use control_data, only:lsecondary
-! use MD_data
-! use MPI_data
-! use map_data
-! use energy
-! use geometry
-! use control
-! 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
-! logical nextp
-!el external nextp
- integer,dimension(isym) :: a
-! parameter(n=symetr)
-!el local variables
- integer :: kkk,i
-
- n=isym
- if (n.eq.1) then
- tabperm(1,1)=1
- return
- endif
- kkk=0
- do i=1,n
- a(i)=i
- enddo
- 10 print *,(a(i),i=1,n)
- kkk=kkk+1
- do i=1,n
- tabperm(kkk,i)=a(i)
-! write (iout,*) "tututu", kkk
- enddo
- if(nextp(n,a)) go to 10
- return
- end subroutine permut
-!-----------------------------------------------------------------------------
- logical function nextp(n,a)
-
- integer :: n,i,j,k,t
-! logical :: nextp
- integer,dimension(n) :: a
- i=n-1
- 10 if(a(i).lt.a(i+1)) go to 20
- i=i-1
- if(i.eq.0) go to 20
- go to 10
- 20 j=i+1
- k=n
- 30 t=a(j)
- a(j)=a(k)
- a(k)=t
- j=j+1
- k=k-1
- if(j.lt.k) go to 30
- j=i
- if(j.ne.0) go to 40
- nextp=.false.
- return
- 40 j=j+1
- if(a(j).lt.a(i)) go to 40
- t=a(i)
- a(i)=a(j)
- a(j)=t
- nextp=.true.
- return
- end function nextp
+#endif
!-----------------------------------------------------------------------------
! printmat.f
!-----------------------------------------------------------------------------
! include 'COMMON.CONTROL'
! include 'COMMON.DISTFIT'
! include 'COMMON.SETUP'
- integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,ity!,&
+ integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift!,ity!,&
! ishift_pdb
logical :: lprn=.true.,fail
real(kind=8),dimension(3) :: e1,e2,e3
allocate(hfrag(2,maxres/3)) !(2,maxres/3)
allocate(bfrag(4,maxres/3)) !(4,maxres/3)
-!elwrite(iout,*)"poczatek read pdb"
-
do i=1,100000
read (ipdbin,'(a80)',end=10) card
! write (iout,'(a)') card
if (card(:5).eq.'HELIX') then
- nhfrag=nhfrag+1
- lsecondary=.true.
- read(card(22:25),*) hfrag(1,nhfrag)
- read(card(34:37),*) hfrag(2,nhfrag)
+ nhfrag=nhfrag+1
+ lsecondary=.true.
+ read(card(22:25),*) hfrag(1,nhfrag)
+ read(card(34:37),*) hfrag(2,nhfrag)
endif
if (card(:5).eq.'SHEET') then
- nbfrag=nbfrag+1
- lsecondary=.true.
- read(card(24:26),*) bfrag(1,nbfrag)
- read(card(35:37),*) bfrag(2,nbfrag)
+ nbfrag=nbfrag+1
+ lsecondary=.true.
+ read(card(24:26),*) bfrag(1,nbfrag)
+ read(card(35:37),*) bfrag(2,nbfrag)
!rc----------------------------------------
!rc to be corrected !!!
- bfrag(3,nbfrag)=bfrag(1,nbfrag)
- bfrag(4,nbfrag)=bfrag(2,nbfrag)
+ bfrag(3,nbfrag)=bfrag(1,nbfrag)
+ bfrag(4,nbfrag)=bfrag(2,nbfrag)
!rc----------------------------------------
endif
if (card(:3).eq.'END') then
goto 10
else if (card(:3).eq.'TER') then
! End current chain
- ires_old=ires+1
+ ires_old=ires+2
ishift1=ishift1+1
itype(ires_old)=ntyp1
+ itype(ires_old-1)=ntyp1
ibeg=2
! write (iout,*) "Chain ended",ires,ishift,ires_old
if (unres_pdb) then
enddo
else
call sccenter(ires,iii,sccor)
+! iii=0
endif
iii=0
endif
! write (iout,*) "Calculating sidechain center iii",iii
if (unres_pdb) then
do j=1,3
- dc(j,ires+nres)=sccor(j,iii)
+ dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
enddo
else
call sccenter(ires_old,iii,sccor)
endif
ires=ires-ishift+ishift1
ires_old=ires
-! write (iout,*) "ishift",ishift," ires",ires,
-! & " ires_old",ires_old
+! write (iout,*) "ishift",ishift," ires",ires,&
+! " ires_old",ires_old
ibeg=0
else if (ibeg.eq.2) then
! Start a new chain
-! ishift=-ires_old+ires-1
-! ishift1=ishift1+1
+ ishift=-ires_old+ires-1 !!!!!
+ ishift1=ishift1-1 !!!!!
! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
ires=ires-ishift+ishift1
ires_old=ires
if (card(27:27).eq."A" .or. card(27:27).eq."B") then
! ishift1=ishift1+1
endif
-! write (2,*) "ires",ires," res ",res," ity",ity
+! write (2,*) "ires",ires," res ",res!," ity"!,ity
if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
res.eq.'NHE'.and.atom(:2).eq.'HN') then
read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
-! write (iout,*) "backbone ",atom
+! write (iout,*) "backbone ",atom
#ifdef DEBUG
write (iout,'(2i3,2x,a,3f8.3)') &
ires,itype(ires),res,(c(j,ires),j=1,3)
nres=ires
do i=2,nres-1
! write (iout,*) i,itype(i)
- if (itype(i).eq.ntyp1) then
+! if (itype(i).eq.ntyp1) then
! write (iout,*) "dummy",i,itype(i)
- do j=1,3
- c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
+! do j=1,3
+! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
! c(j,i)=(c(j,i-1)+c(j,i+1))/2
- dc(j,i)=c(j,i)
- enddo
- endif
+! dc(j,i)=c(j,i)
+! enddo
+! endif
+ if (itype(i).eq.ntyp1) then
+ if (itype(i+1).eq.ntyp1) then
+! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
+! first is connected prevous chain (itype(i+1).eq.ntyp1)=true
+! second dummy atom is conected to next chain itype(i+1).eq.ntyp1=false
+ if (unres_pdb) then
+! 2/15/2013 by Adam: corrected insertion of the last dummy residue
+! print *,i,'tu dochodze'
+ call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif !fail
+ print *,i,'a tu?'
+ do j=1,3
+ c(j,i)=c(j,i-1)-1.9d0*e2(j)
+ enddo
+ else !unres_pdb
+ do j=1,3
+ dcj=(c(j,i-2)-c(j,i-3))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,i)=c(j,i-1)+dcj
+ c(j,nres+i)=c(j,i)
+ enddo
+ endif !unres_pdb
+ else !itype(i+1).eq.ntyp1
+ if (unres_pdb) then
+! 2/15/2013 by Adam: corrected insertion of the first dummy residue
+ call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
+ if (fail) then
+ e2(1)=0.0d0
+ e2(2)=1.0d0
+ e2(3)=0.0d0
+ endif
+ do j=1,3
+ c(j,i)=c(j,i+1)-1.9d0*e2(j)
+ enddo
+ else !unres_pdb
+ do j=1,3
+ dcj=(c(j,i+3)-c(j,i+2))/2.0
+ if (dcj.eq.0) dcj=1.23591524223
+ c(j,i)=c(j,i+1)-dcj
+ c(j,nres+i)=c(j,i)
+ enddo
+ endif !unres_pdb
+ endif !itype(i+1).eq.ntyp1
+ endif !itype.eq.ntyp1
+
enddo
! Calculate the CM of the last side chain.
if (iii.gt.0) then
e2(3)=0.0d0
endif
do j=1,3
- c(j,nres)=c(j,nres-1)-3.8d0*e2(j)
+ c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
enddo
else
do j=1,3
- dcj=c(j,nres-2)-c(j,nres-3)
+ dcj=(c(j,nres-2)-c(j,nres-3))/2.0
c(j,nres)=c(j,nres-1)+dcj
c(j,2*nres)=c(j,nres)
enddo
endif
endif
+!el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
+#ifdef WHAM_RUN
+ if (nres.ne.nres0) then
+ write (iout,*) "Error: wrong parameter value: NRES=",nres,&
+ " NRES0=",nres0
+ stop "Error nres value in WHAM input"
+ endif
+#endif
!---------------------------------
!el reallocate tables
- do i=1,maxres/3
- do j=1,2
- hfrag_alloc(j,i)=hfrag(j,i)
- enddo
- do j=1,4
- bfrag_alloc(j,i)=bfrag(j,i)
- enddo
- enddo
+! do i=1,maxres/3
+! do j=1,2
+! hfrag_alloc(j,i)=hfrag(j,i)
+! enddo
+! do j=1,4
+! bfrag_alloc(j,i)=bfrag(j,i)
+! enddo
+! enddo
- deallocate(hfrag)
- deallocate(bfrag)
- allocate(hfrag(2,nres/3)) !(2,maxres/3)
+! deallocate(hfrag)
+! deallocate(bfrag)
+! allocate(hfrag(2,nres/3)) !(2,maxres/3)
!el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
!el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
- allocate(bfrag(4,nres/3)) !(4,maxres/3)
+! allocate(bfrag(4,nres/3)) !(4,maxres/3)
- do i=1,nhfrag
- do j=1,2
- hfrag(j,i)=hfrag_alloc(j,i)
- enddo
- enddo
- do i=1,nbfrag
- do j=1,4
- bfrag(j,i)=bfrag_alloc(j,i)
- enddo
- enddo
+! do i=1,nhfrag
+! do j=1,2
+! hfrag(j,i)=hfrag_alloc(j,i)
+! enddo
+! enddo
+! do i=1,nbfrag
+! do j=1,4
+! bfrag(j,i)=bfrag_alloc(j,i)
+! enddo
+! enddo
!el end reallocate tables
!---------------------------------
do i=2,nres-1
e2(3)=0.0d0
endif
do j=1,3
- c(j,1)=c(j,2)-3.8d0*e2(j)
+ c(j,1)=c(j,2)-1.9d0*e2(j)
enddo
else
do j=1,3
- dcj=c(j,4)-c(j,3)
+ dcj=(c(j,4)-c(j,3))/2.0
c(j,1)=c(j,2)-dcj
c(j,nres+1)=c(j,1)
enddo
(c(j,ires+nres),j=1,3)
enddo
endif
+! znamy już nres wiec mozna alokowac tablice
! Calculate internal coordinates.
if(me.eq.king.or..not.out1file)then
write (iout,'(a)') &
enddo
endif
- if(.not.allocated(vbld)) allocate(vbld(2*nres))
- if(.not.allocated(vbld_inv)) allocate(vbld_inv(2*nres))
- if(.not.allocated(theta)) allocate(theta(nres+2))
+ if(.not.allocated(vbld)) then
+ allocate(vbld(2*nres))
+ do i=1,2*nres
+ vbld(i)=0.d0
+ enddo
+ endif
+ if(.not.allocated(vbld_inv)) then
+ allocate(vbld_inv(2*nres))
+ do i=1,2*nres
+ vbld_inv(i)=0.d0
+ enddo
+ endif
+!!!el
+ if(.not.allocated(theta)) then
+ allocate(theta(nres+2))
+ theta(:)=0.0d0
+ endif
+
if(.not.allocated(phi)) allocate(phi(nres+2))
if(.not.allocated(alph)) allocate(alph(nres+2))
if(.not.allocated(omeg)) allocate(omeg(nres+2))
- if(.not.allocated(theta)) allocate(theta(nres+2))
if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
if(.not.allocated(phiref)) allocate(phiref(nres+2))
if(.not.allocated(costtab)) allocate(costtab(nres))
if(.not.allocated(xxref)) allocate(xxref(nres))
if(.not.allocated(yyref)) allocate(yyref(nres))
if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
- if(.not.allocated(theta)) allocate(theta(nres))
- if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres))
- if(.not.allocated(theta)) allocate(theta(nres))
+ if(.not.allocated(dc_norm)) then
+! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
+ allocate(dc_norm(3,0:2*nres+2))
+ dc_norm(:,:)=0.d0
+ endif
+
call int_from_cart(.true.,.false.)
- call sc_loc_geom(.true.)
-! wczesbiej bylo false
+ call sc_loc_geom(.false.)
do i=1,nres
thetaref(i)=theta(i)
phiref(i)=phi(i)
enddo
+! do i=1,2*nres
+! vbld_inv(i)=0.d0
+! vbld(i)=0.d0
+! enddo
+
do i=1,nres-1
do j=1,3
dc(j,i)=c(j,i+1)-c(j,i)
dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
enddo
enddo
-
do i=2,nres-1
do j=1,3
dc(j,i+nres)=c(j,i+nres)-c(j,i)
! enddo
! enddo
!
- allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
- allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
- allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
+ if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
+ if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
+ if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
!-----------------------------
kkk=1
lll=0
! enddiagnostic
! makes copy of chains
write (iout,*) "symetr", symetr
+ do j=1,3
+ dc(j,0)=c(j,1)
+ enddo
if (symetr.gt.1) then
call permut(symetr)
enddo
enddo
ishift_pdb=ishift
-!---------------------
-! el reallocate array
- do i=1,2*nres+2
- do kkk=1,nperm
- cref_alloc(1,i,kkk)=cref(1,i,kkk)
- cref_alloc(2,i,kkk)=cref(2,i,kkk)
- cref_alloc(3,i,kkk)=cref(3,i,kkk)
- enddo
- enddo
-!el deallocate(cref)
-!el allocate(cref(3,2*nres+2,nperm)) !(3,maxres2+2,maxperm)
-
- do i=1,2*nres+2
- do kkk=1,nperm
- cref(1,i,kkk)=cref_alloc(1,i,kkk)
- cref(2,i,kkk)=cref_alloc(2,i,kkk)
- cref(3,i,kkk)=cref_alloc(3,i,kkk)
- enddo
- enddo
-!---------------------
+
return
end subroutine readpdb
-#ifndef WHAM_RUN
+#if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
!-----------------------------------------------------------------------------
! readrtns_CSA.F
!-----------------------------------------------------------------------------
character(len=640) :: controlcard
real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
-
+ integer i
nglob_csa=0
eglob_csa=1d99
timem=timlim
modecalc=0
call reada(controlcard,"T_BATH",t_bath,300.0d0)
+!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
+!C solvent accesible area
+!C 1 the newly introduced function
+!C 2 reseved for further possible developement
+ call readi(controlcard,'SHIELD',shield_mode,0)
+!C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write(iout,*) "shield_mode",shield_mode
+!C Varibles set size of box
+ call reada(controlcard,'BOXX',boxxsize,100.0d0)
+ call reada(controlcard,'BOXY',boxysize,100.0d0)
+ call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+! CUTOFFF ON ELECTROSTATICS
+ call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
+ call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
+ write(iout,*) "R_CUT_ELE=",r_cut_ele
+! Lipidic parameters
+ call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
+ call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
+ if (lipthick.gt.0.0d0) then
+ bordliptop=(boxzsize+lipthick)/2.0
+ bordlipbot=bordliptop-lipthick
+ if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
+ write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
+ buflipbot=bordlipbot+lipbufthick
+ bufliptop=bordliptop-lipbufthick
+ if ((lipbufthick*2.0d0).gt.lipthick) &
+ write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
+ endif !lipthick.gt.0
+ write(iout,*) "bordliptop=",bordliptop
+ write(iout,*) "bordlipbot=",bordlipbot
+ write(iout,*) "bufliptop=",bufliptop
+ write(iout,*) "buflipbot=",buflipbot
+ write (iout,*) "SHIELD MODE",shield_mode
+
+!C-------------------------
minim=(index(controlcard,'MINIMIZE').gt.0)
dccart=(index(controlcard,'CART').gt.0)
overlapsc=(index(controlcard,'OVERLAP').gt.0)
if(me.eq.king.or..not.out1file) &
write (iout,'(2a)') diagmeth(kdiag),&
' routine used to diagonalize matrices.'
+ if (shield_mode.gt.0) then
+ pi=3.141592d0
+!C VSolvSphere the volume of solving sphere
+!C print *,pi,"pi"
+!C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
+!C there will be no distinction between proline peptide group and normal peptide
+!C group in case of shielding parameters
+ VSolvSphere=4.0/3.0*pi*rpp(1,1)**3
+ VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(rpp(1,1)/2.0)**3
+ write (iout,*) VSolvSphere,VSolvSphere_div
+!C long axis of side chain
+ do i=1,ntyp
+ long_r_sidechain(i)=vbldsc0(1,i)
+ short_r_sidechain(i)=sigma0(i)
+ write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
+ sigma0(i)
+ enddo
+ buff_shield=1.0d0
+ endif
return
end subroutine read_control
!-----------------------------------------------------------------------------
!-----------------------------------------------------------------------------
subroutine openunits
- use energy_data, only: usampl
+ use MD_data, only: usampl
use csa_data
use MPI_data
-! use MD
use control_data, only:out1file
use control, only: getenv_loc
! implicit real*8 (a-h,o-z)
! print *,"Processor",myrank," opened file IELEP"
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old',action='read')
+ call getenv_loc('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old',action='read')
! print *,"Processor",myrank," opened file ISIDEP"
! print *,"Processor",myrank," opened parameter files"
#elif (defined G77)
open (ielep,file=elename,status='old')
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old')
+ call getenv_loc('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old')
#else
open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
readonly)
open (ielep,file=elename,status='old',readonly)
call getenv_loc('SIDEPAR',sidename)
open (isidep,file=sidename,status='old',readonly)
+ call getenv_loc('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old',action='read')
#ifndef CRYST_SC
call getenv_loc('ROTPARPDB',rotname_pdb)
open (irotam_pdb,file=rotname_pdb,status='old',action='read')
thetname(:ilen(thetname))
write (iout,*) "Rotamer parameter file : ",&
rotname(:ilen(rotname))
+!el----
+#ifndef CRYST_THETA
write (iout,*) "Thetpdb parameter file : ",&
thetname_pdb(:ilen(thetname_pdb))
+#endif
+!el
write (iout,*) "Threading database : ",&
patname(:ilen(patname))
if (lentmp.ne.0) &
subroutine readrst
use geometry_data, only: nres,dc
- use energy_data, only: usampl,iset
use MD_data
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
open(irest2,file=rest2name,status='unknown')
read(irest2,*) totT,EK,potE,totE,t_bath
- do i=1,2*nres
+! do i=1,2*nres
+! AL 4/17/17: Now reading d_t(0,:) too
+ do i=0,2*nres
read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
enddo
- do i=1,2*nres
+! do i=1,2*nres
+! AL 4/17/17: Now reading d_c(0,:) too
+ do i=0,2*nres
read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
enddo
if(usampl) then
projects / unres4.git / blobdiff