! Get parameter filenames and open the parameter files.
call mygetenv('BONDPAR',bondname)
open (ibond,file=bondname,status='old')
+ call mygetenv('BONDPAR_NUCL',bondname_nucl)
+ open (ibond_nucl,file=bondname_nucl,status='old')
call mygetenv('THETPAR',thetname)
open (ithep,file=thetname,status='old')
call mygetenv('ROTPAR',rotname)
open (isidep,file=sidename,status='old')
call mygetenv('SIDEP',sidepname)
open (isidep1,file=sidepname,status="old")
+ call mygetenv('THETPAR_NUCL',thetname_nucl)
+ open (ithep_nucl,file=thetname_nucl,status='old')
+ call mygetenv('ROTPAR_NUCL',rotname_nucl)
+ open (irotam_nucl,file=rotname_nucl,status='old')
+ call mygetenv('TORPAR_NUCL',torname_nucl)
+ open (itorp_nucl,file=torname_nucl,status='old')
+ call mygetenv('TORDPAR_NUCL',tordname_nucl)
+ open (itordp_nucl,file=tordname_nucl,status='old')
+ call mygetenv('SIDEPAR_NUCL',sidename_nucl)
+ open (isidep_nucl,file=sidename_nucl,status='old')
+ call mygetenv('SIDEPAR_SCBASE',sidename_scbase)
+ open (isidep_scbase,file=sidename_scbase,status='old')
+ call mygetenv('PEPPAR_PEPBASE',pepname_pepbase)
+ open (isidep_pepbase,file=pepname_pepbase,status='old')
+ call mygetenv('SCPAR_PHOSPH',pepname_scpho)
+ open (isidep_scpho,file=pepname_scpho,status='old')
+ call mygetenv('PEPPAR_PHOSPH',pepname_peppho)
+ open (isidep_peppho,file=pepname_peppho,status='old')
+ call mygetenv('LIPTRANPAR',liptranname)
+ open (iliptranpar,file=liptranname,status='old',action='read')
+ call mygetenv('TUBEPAR',tubename)
+ open (itube,file=tubename,status='old',action='read')
+ call mygetenv('IONPAR',ionname)
+ open (iion,file=ionname,status='old',action='read')
+
+ call mygetenv('SCPPAR_NUCL',scpname_nucl)
+ open (iscpp_nucl,file=scpname_nucl,status='old')
+
+
#ifndef OLDSCP
!
! 8/9/01 In the newest version SCp interaction constants are read from a file
use geometry_data
use energy_data
use control_data, only: maxterm,maxlor,maxterm_sccor,& !maxtor
- maxtermd_1,maxtermd_2 !,maxthetyp,maxthetyp1
+ maxtermd_1,maxtermd_2,tor_mode,scelemode !,maxthetyp,maxthetyp1
use MD_data
!el use MPI_data
!el use map_data
character(len=1) :: t1,t2,t3
character(len=1) :: onelett(4) = (/"G","A","P","D"/)
character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
- logical :: lprint
+ logical :: lprint,SPLIT_FOURIERTOR
real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
character(len=800) :: controlcard
character(len=256) :: bondname_t,thetname_t,rotname_t,torname_t,&
!el integer ilen
!el external ilen
character(len=16) :: key
- integer :: iparm
+ integer :: iparm,nkcctyp
!el real(kind=8) :: ip,mp
real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,epsijlip,&
- sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm
+ sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm, &
+ epspeptube,epssctube,sigmapeptube, &
+ sigmasctube
+
real(kind=8) :: v0ij,v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,rjunk,&
res1
- integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n
+ integer :: i,j,ichir1,ichir2,k,l,m,kk,ii,mm,junk,lll,ll,llll,n,jj
integer :: nlobi,iblock,maxinter,iscprol,ncatprotparm
+ character*3 string
+
!
! Body
!
wscloc=ww(12)
wtor=ww(13)
wtor_d=ww(14)
+ wstrain=ww(15)
wvdwpp=ww(16)
wbond=ww(18)
wsccor=ww(19)
- wcatcat=ww(44)
- wcatprot=ww(43)
+ wcatcat=ww(42)
+ wcatprot=ww(41)
+ wcorr3_nucl=ww(38)
+ wcorr_nucl=ww(37)
+ wtor_d_nucl=ww(36)
+ wtor_nucl=ww(35)
+ wsbloc=ww(34)
+ wang_nucl=ww(33)
+ wbond_nucl=ww(32)
+ welsb=ww(31)
+ wvdwsb=ww(30)
+ welpsb=ww(29)
+ wvdwpsb=ww(28)
+ welpp=ww(27)
+ wvdwpp_nucl=ww(26)
+ wscbase=ww(46)
+ wpepbase=ww(47)
+ wscpho=ww(48)
+ wpeppho=ww(49)
+! print *,"KURWA",ww(48)
+! "WSCBASE ","WPEPBASE ","WSCPHO ","WPEPPHO "
+! "WVDWPP ","WELPP ","WVDWPSB ","WELPSB ","WVDWSB ",&
+! "WELSB ","WBOND ","WANG ","WSBLOC ","WTOR ",&
+! "WTORD ","WCORR ","WCORR3 ","WNULL ","WNULL ",&
+! "WCATPROT ","WCATCAT
+! +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+! +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+! +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+! +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
endif
!
weights(12)=wscloc
weights(13)=wtor
weights(14)=wtor_d
- weights(15)=0 !wstrain !
- weights(16)=0 !wvdwpp !
+ weights(15)=wstrain !0
+ weights(16)=wvdwpp !
weights(17)=wbond
weights(18)=0 !scal14 !
weights(21)=wsccor
- weights(43)=wcatprot
- weights(44)=wcatcat
-
+ weights(42)=wcatprot
+ weights(41)=wcatcat
+ weights(26)= wvdwpp_nucl
+
+ weights(27) =welpp
+ weights(28) =wvdwpsb
+ weights(29) =welpsb
+ weights(30) =wvdwsb
+ weights(31) =welsb
+ weights(32) =wbond_nucl
+ weights(33) =wang_nucl
+ weights(34) =wsbloc
+ weights(35) =wtor_nucl
+ weights(36) =wtor_d_nucl
+ weights(37) =wcorr_nucl
+ weights(38) =wcorr3_nucl
+ weights(41) =wcatcat
+ weights(42) =wcatprot
+ weights(46) =wscbase
+ weights(47) =wscpho
+ weights(48) =wpeppho
! el--------
call card_concat(controlcard,.false.)
allocate(nbondterm(ntyp)) !(ntyp)
allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
+ allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
+ allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
+ allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
+
!el allocate(msc(ntyp+1)) !(ntyp+1)
!el allocate(isc(ntyp+1)) !(ntyp+1)
!el allocate(restok(ntyp+1)) !(ntyp+1)
read (iion,*) (catprm(i,k),i=1,ncatprotparm)
enddo
print *, catprm
+ read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
+ do i=1,ntyp_molec(2)
+ nbondterm_nucl(i)=1
+ read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),rjunk,restok(i,2)
+! dsc(i) = vbldsc0_nucl(1,i)
+! if (i.eq.10) then
+! dsc_inv(i)=0.0D0
+! else
+! dsc_inv(i)=1.0D0/dsc(i)
+! endif
+ enddo
!----------------------------------------------------
! Read the parameters of Utheta determined from ab initio surfaces
! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
!
+ allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
! write (iout,*) "tu dochodze"
+ IF (tor_mode.eq.0) THEN
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
read (ithep,*) nthetyp,ntheterm,ntheterm2,&
ntheterm3,nsingle,ndouble
nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
!----------------------------------------------------
- allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
+! allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
allocate(aa0thet(-nthetyp-1:nthetyp+1,&
-nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
!(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
!
! Control printout of the coefficients of virtual-bond-angle potentials
!
-do iblock=1,2
+ do iblock=1,2
if (lprint) then
write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
do i=1,nthetyp+1
enddo
call flush(iout)
endif
-enddo
+ enddo
+ ELSE
+!C here will be the apropriate recalibrating for D-aminoacid
+ read (ithep,*) nthetyp
+ allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
+ allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
+ do i=-nthetyp+1,nthetyp-1
+ read (ithep,*) nbend_kcc_Tb(i)
+ do j=0,nbend_kcc_Tb(i)
+ read (ithep,*) ijunk,v1bend_chyb(j,i)
+ enddo
+ enddo
+ if (lprint) then
+ write (iout,'(a)') &
+ "Parameters of the valence-only potentials"
+ do i=-nthetyp+1,nthetyp-1
+ write (iout,'(2a)') "Type ",toronelet(i)
+ do k=0,nbend_kcc_Tb(i)
+ write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
+ enddo
+ enddo
+ endif
+ ENDIF ! TOR_MODE
+
#endif
+!--------------- Reading theta parameters for nucleic acid-------
+ read (ithep_nucl,*) nthetyp_nucl,ntheterm_nucl,&
+ ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
+ nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
+ allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
+ allocate(aa0thet_nucl(nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+!(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+ allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+!(maxtheterm,-maxthetyp1:maxthetyp1,&
+! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+ allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+ allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+ allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+ allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1))
+!(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
+! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
+ allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
+ allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
+ nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
+
+!(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
+! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
+
+ read (ithep_nucl,*) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
+
+ aa0thet_nucl(:,:,:)=0.0d0
+ aathet_nucl(:,:,:,:)=0.0d0
+ bbthet_nucl(:,:,:,:,:)=0.0d0
+ ccthet_nucl(:,:,:,:,:)=0.0d0
+ ddthet_nucl(:,:,:,:,:)=0.0d0
+ eethet_nucl(:,:,:,:,:)=0.0d0
+ ffthet_nucl(:,:,:,:,:,:)=0.0d0
+ ggthet_nucl(:,:,:,:,:,:)=0.0d0
+
+ do i=1,nthetyp_nucl
+ do j=1,nthetyp_nucl
+ do k=1,nthetyp_nucl
+ read (ithep_nucl,'(3a)') t1,t2,t3
+ read (ithep_nucl,*) aa0thet_nucl(i,j,k)
+ read (ithep_nucl,*)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
+ read (ithep_nucl,*) &
+ (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
+ (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
+ (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
+ (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
+ read (ithep_nucl,*) &
+ (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
+ ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
+ llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
+ enddo
+ enddo
+ enddo
+
+
!-------------------------------------------
allocate(nlob(ntyp1)) !(ntyp1)
allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
enddo
#endif
close(irotam)
+!---------reading nucleic acid parameters for rotamers-------------------
+ allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
+ do i=1,ntyp_molec(2)
+ read (irotam_nucl,*)
+ do j=1,9
+ read(irotam_nucl,*) sc_parmin_nucl(j,i)
+ enddo
+ enddo
+ close(irotam_nucl)
+ if (lprint) then
+ write (iout,*)
+ write (iout,*) "Base rotamer parameters"
+ do i=1,ntyp_molec(2)
+ write (iout,'(a)') restyp(i,2)
+ write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
+ enddo
+ endif
+
+
+ read (ifourier,*) nloctyp
+!el write(iout,*)"nloctyp",nloctyp
+ SPLIT_FOURIERTOR = nloctyp.lt.0
+ nloctyp = iabs(nloctyp)
+#ifdef NEWCORR
+ if (.not.allocated(itype2loc)) allocate(itype2loc(-ntyp1:ntyp1))
+ print *,"shape",shape(itype2loc)
+ allocate(iloctyp(-nloctyp:nloctyp))
+ allocate(bnew1(3,2,-nloctyp:nloctyp))
+ allocate(bnew2(3,2,-nloctyp:nloctyp))
+ allocate(ccnew(3,2,-nloctyp:nloctyp))
+ allocate(ddnew(3,2,-nloctyp:nloctyp))
+ allocate(e0new(3,-nloctyp:nloctyp))
+ allocate(eenew(2,2,2,-nloctyp:nloctyp))
+ allocate(bnew1tor(3,2,-nloctyp:nloctyp))
+ allocate(bnew2tor(3,2,-nloctyp:nloctyp))
+ allocate(ccnewtor(3,2,-nloctyp:nloctyp))
+ allocate(ddnewtor(3,2,-nloctyp:nloctyp))
+ allocate(e0newtor(3,-nloctyp:nloctyp))
+ allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
+
+ read (ifourier,*) (itype2loc(i),i=1,ntyp)
+ read (ifourier,*) (iloctyp(i),i=0,nloctyp-1)
+ itype2loc(ntyp1)=nloctyp
+ iloctyp(nloctyp)=ntyp1
+ do i=1,ntyp1
+ itype2loc(-i)=-itype2loc(i)
+ enddo
+#else
+ allocate(iloctyp(-nloctyp:nloctyp))
+ allocate(itype2loc(-ntyp1:ntyp1))
+ iloctyp(0)=10
+ iloctyp(1)=9
+ iloctyp(2)=20
+ iloctyp(3)=ntyp1
+#endif
+ do i=1,nloctyp
+ iloctyp(-i)=-iloctyp(i)
+ enddo
+!c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
+!c write (iout,*) "nloctyp",nloctyp,
+!c & " iloctyp",(iloctyp(i),i=0,nloctyp)
+!c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
+!c write (iout,*) "nloctyp",nloctyp,
+!c & " iloctyp",(iloctyp(i),i=0,nloctyp)
+#ifdef NEWCORR
+ do i=0,nloctyp-1
+!c write (iout,*) "NEWCORR",i
+ read (ifourier,*)
+ do ii=1,3
+ do j=1,2
+ read (ifourier,*) bnew1(ii,j,i)
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR BNEW1"
+!c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
+ do ii=1,3
+ do j=1,2
+ read (ifourier,*) bnew2(ii,j,i)
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR BNEW2"
+!c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
+ do kk=1,3
+ read (ifourier,*) ccnew(kk,1,i)
+ read (ifourier,*) ccnew(kk,2,i)
+ enddo
+!c write (iout,*) "NEWCORR CCNEW"
+!c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
+ do kk=1,3
+ read (ifourier,*) ddnew(kk,1,i)
+ read (ifourier,*) ddnew(kk,2,i)
+ enddo
+!c write (iout,*) "NEWCORR DDNEW"
+!c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
+ do ii=1,2
+ do jj=1,2
+ do kk=1,2
+ read (ifourier,*) eenew(ii,jj,kk,i)
+ enddo
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR EENEW1"
+!c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
+ do ii=1,3
+ read (ifourier,*) e0new(ii,i)
+ enddo
+!c write (iout,*) (e0new(ii,i),ii=1,3)
+ enddo
+!c write (iout,*) "NEWCORR EENEW"
+ do i=0,nloctyp-1
+ do ii=1,3
+ ccnew(ii,1,i)=ccnew(ii,1,i)/2
+ ccnew(ii,2,i)=ccnew(ii,2,i)/2
+ ddnew(ii,1,i)=ddnew(ii,1,i)/2
+ ddnew(ii,2,i)=ddnew(ii,2,i)/2
+ enddo
+ enddo
+ do i=1,nloctyp-1
+ do ii=1,3
+ bnew1(ii,1,-i)= bnew1(ii,1,i)
+ bnew1(ii,2,-i)=-bnew1(ii,2,i)
+ bnew2(ii,1,-i)= bnew2(ii,1,i)
+ bnew2(ii,2,-i)=-bnew2(ii,2,i)
+ enddo
+ do ii=1,3
+!c ccnew(ii,1,i)=ccnew(ii,1,i)/2
+!c ccnew(ii,2,i)=ccnew(ii,2,i)/2
+!c ddnew(ii,1,i)=ddnew(ii,1,i)/2
+!c ddnew(ii,2,i)=ddnew(ii,2,i)/2
+ ccnew(ii,1,-i)=ccnew(ii,1,i)
+ ccnew(ii,2,-i)=-ccnew(ii,2,i)
+ ddnew(ii,1,-i)=ddnew(ii,1,i)
+ ddnew(ii,2,-i)=-ddnew(ii,2,i)
+ enddo
+ e0new(1,-i)= e0new(1,i)
+ e0new(2,-i)=-e0new(2,i)
+ e0new(3,-i)=-e0new(3,i)
+ do kk=1,2
+ eenew(kk,1,1,-i)= eenew(kk,1,1,i)
+ eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
+ eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
+ eenew(kk,2,2,-i)= eenew(kk,2,2,i)
+ enddo
+ enddo
+ if (lprint) then
+ write (iout,'(a)') "Coefficients of the multibody terms"
+ do i=-nloctyp+1,nloctyp-1
+ write (iout,*) "Type: ",onelet(iloctyp(i))
+ write (iout,*) "Coefficients of the expansion of B1"
+ do j=1,2
+ write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
+ enddo
+ write (iout,*) "Coefficients of the expansion of B2"
+ do j=1,2
+ write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
+ enddo
+ write (iout,*) "Coefficients of the expansion of C"
+ write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
+ write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
+ write (iout,*) "Coefficients of the expansion of D"
+ write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
+ write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
+ write (iout,*) "Coefficients of the expansion of E"
+ write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
+ do j=1,2
+ do k=1,2
+ write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
+ enddo
+ enddo
+ enddo
+ endif
+ IF (SPLIT_FOURIERTOR) THEN
+ do i=0,nloctyp-1
+!c write (iout,*) "NEWCORR TOR",i
+ read (ifourier,*)
+ do ii=1,3
+ do j=1,2
+ read (ifourier,*) bnew1tor(ii,j,i)
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR BNEW1 TOR"
+!c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
+ do ii=1,3
+ do j=1,2
+ read (ifourier,*) bnew2tor(ii,j,i)
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR BNEW2 TOR"
+!c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
+ do kk=1,3
+ read (ifourier,*) ccnewtor(kk,1,i)
+ read (ifourier,*) ccnewtor(kk,2,i)
+ enddo
+!c write (iout,*) "NEWCORR CCNEW TOR"
+!c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
+ do kk=1,3
+ read (ifourier,*) ddnewtor(kk,1,i)
+ read (ifourier,*) ddnewtor(kk,2,i)
+ enddo
+!c write (iout,*) "NEWCORR DDNEW TOR"
+!c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
+ do ii=1,2
+ do jj=1,2
+ do kk=1,2
+ read (ifourier,*) eenewtor(ii,jj,kk,i)
+ enddo
+ enddo
+ enddo
+!c write (iout,*) "NEWCORR EENEW1 TOR"
+!c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
+ do ii=1,3
+ read (ifourier,*) e0newtor(ii,i)
+ enddo
+!c write (iout,*) (e0newtor(ii,i),ii=1,3)
+ enddo
+!c write (iout,*) "NEWCORR EENEW TOR"
+ do i=0,nloctyp-1
+ do ii=1,3
+ ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
+ ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
+ ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
+ ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
+ enddo
+ enddo
+ do i=1,nloctyp-1
+ do ii=1,3
+ bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
+ bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
+ bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
+ bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
+ enddo
+ do ii=1,3
+ ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
+ ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
+ ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
+ ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
+ enddo
+ e0newtor(1,-i)= e0newtor(1,i)
+ e0newtor(2,-i)=-e0newtor(2,i)
+ e0newtor(3,-i)=-e0newtor(3,i)
+ do kk=1,2
+ eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
+ eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
+ eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
+ eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
+ enddo
+ enddo
+ if (lprint) then
+ write (iout,'(a)') &
+ "Single-body coefficients of the torsional potentials"
+ do i=-nloctyp+1,nloctyp-1
+ write (iout,*) "Type: ",onelet(iloctyp(i))
+ write (iout,*) "Coefficients of the expansion of B1tor"
+ do j=1,2
+ write (iout,'(3hB1(,i1,1h),3f10.5)') &
+ j,(bnew1tor(k,j,i),k=1,3)
+ enddo
+ write (iout,*) "Coefficients of the expansion of B2tor"
+ do j=1,2
+ write (iout,'(3hB2(,i1,1h),3f10.5)') &
+ j,(bnew2tor(k,j,i),k=1,3)
+ enddo
+ write (iout,*) "Coefficients of the expansion of Ctor"
+ write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
+ write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
+ write (iout,*) "Coefficients of the expansion of Dtor"
+ write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
+ write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
+ write (iout,*) "Coefficients of the expansion of Etor"
+ write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
+ do j=1,2
+ do k=1,2
+ write (iout,'(1hE,2i1,2f10.5)') &
+ j,k,(eenewtor(l,j,k,i),l=1,2)
+ enddo
+ enddo
+ enddo
+ endif
+ ELSE
+ do i=-nloctyp+1,nloctyp-1
+ do ii=1,3
+ do j=1,2
+ bnew1tor(ii,j,i)=bnew1(ii,j,i)
+ enddo
+ enddo
+ do ii=1,3
+ do j=1,2
+ bnew2tor(ii,j,i)=bnew2(ii,j,i)
+ enddo
+ enddo
+ do ii=1,3
+ ccnewtor(ii,1,i)=ccnew(ii,1,i)
+ ccnewtor(ii,2,i)=ccnew(ii,2,i)
+ ddnewtor(ii,1,i)=ddnew(ii,1,i)
+ ddnewtor(ii,2,i)=ddnew(ii,2,i)
+ enddo
+ enddo
+ ENDIF !SPLIT_FOURIER_TOR
+#else
+ allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
+ allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
+ allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
+ allocate(b(13,-nloctyp-1:nloctyp+1))
+ if (lprint) &
+ write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
+ do i=0,nloctyp-1
+ read (ifourier,*)
+ read (ifourier,*) (b(ii,i),ii=1,13)
+ if (lprint) then
+ write (iout,*) 'Type ',onelet(iloctyp(i))
+ write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
+ endif
+ if (i.gt.0) then
+ b(2,-i)= b(2,i)
+ b(3,-i)= b(3,i)
+ b(4,-i)=-b(4,i)
+ b(5,-i)=-b(5,i)
+ endif
+ CCold(1,1,i)= b(7,i)
+ CCold(2,2,i)=-b(7,i)
+ CCold(2,1,i)= b(9,i)
+ CCold(1,2,i)= b(9,i)
+ CCold(1,1,-i)= b(7,i)
+ CCold(2,2,-i)=-b(7,i)
+ CCold(2,1,-i)=-b(9,i)
+ CCold(1,2,-i)=-b(9,i)
+ DDold(1,1,i)= b(6,i)
+ DDold(2,2,i)=-b(6,i)
+ DDold(2,1,i)= b(8,i)
+ DDold(1,2,i)= b(8,i)
+ DDold(1,1,-i)= b(6,i)
+ DDold(2,2,-i)=-b(6,i)
+ DDold(2,1,-i)=-b(8,i)
+ DDold(1,2,-i)=-b(8,i)
+ EEold(1,1,i)= b(10,i)+b(11,i)
+ EEold(2,2,i)=-b(10,i)+b(11,i)
+ EEold(2,1,i)= b(12,i)-b(13,i)
+ EEold(1,2,i)= b(12,i)+b(13,i)
+ EEold(1,1,-i)= b(10,i)+b(11,i)
+ EEold(2,2,-i)=-b(10,i)+b(11,i)
+ EEold(2,1,-i)=-b(12,i)+b(13,i)
+ EEold(1,2,-i)=-b(12,i)-b(13,i)
+ write(iout,*) "TU DOCHODZE"
+ print *,"JESTEM"
+ enddo
+ if (lprint) then
+ write (iout,*)
+ write (iout,*) &
+ "Coefficients of the cumulants (independent of valence angles)"
+ do i=-nloctyp+1,nloctyp-1
+ write (iout,*) 'Type ',onelet(iloctyp(i))
+ write (iout,*) 'B1'
+ write(iout,'(2f10.5)') B(3,i),B(5,i)
+ write (iout,*) 'B2'
+ write(iout,'(2f10.5)') B(2,i),B(4,i)
+ write (iout,*) 'CC'
+ do j=1,2
+ write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
+ enddo
+ write(iout,*) 'DD'
+ do j=1,2
+ write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
+ enddo
+ write(iout,*) 'EE'
+ do j=1,2
+ write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
+ enddo
+ enddo
+ endif
+#endif
#ifdef CRYST_TOR
!
! Read torsional parameters in old format
!
! Read torsional parameters
!
+ IF (TOR_MODE.eq.0) THEN
+
allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
read (itorp,*) ntortyp
enddo
enddo
endif
+#ifndef NEWCORR
+ do i=1,ntyp1
+ itype2loc(i)=itortyp(i)
+ enddo
+#endif
+ ELSE IF (TOR_MODE.eq.1) THEN
+
+!C read valence-torsional parameters
+ read (itorp,*) ntortyp
+ nkcctyp=ntortyp
+ write (iout,*) "Valence-torsional parameters read in ntortyp",&
+ ntortyp
+ read (itorp,*) (itortyp(i),i=1,ntyp)
+ write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
+#ifndef NEWCORR
+ do i=1,ntyp1
+ itype2loc(i)=itortyp(i)
+ enddo
+#endif
+ do i=-ntyp,-1
+ itortyp(i)=-itortyp(-i)
+ enddo
+ do i=-ntortyp+1,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
+!C first we read the cos and sin gamma parameters
+ read (itorp,'(13x,a)') string
+ write (iout,*) i,j,string
+ read (itorp,*) &
+ nterm_kcc(j,i),nterm_kcc_Tb(j,i)
+!C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
+ do k=1,nterm_kcc(j,i)
+ do l=1,nterm_kcc_Tb(j,i)
+ do ll=1,nterm_kcc_Tb(j,i)
+ read (itorp,*) ii,jj,kk, &
+ v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ ELSE
+#ifdef NEWCORR
+!c AL 4/8/16: Calculate coefficients from one-body parameters
+ ntortyp=nloctyp
+ allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
+ allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
+ allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
+ allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
+ allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
+
+ do i=-ntyp1,ntyp1
+ print *,i,itortyp(i)
+ itortyp(i)=itype2loc(i)
+ enddo
+ write (iout,*) &
+ "Val-tor parameters calculated from cumulant coefficients ntortyp"&
+ ,ntortyp
+ do i=-ntortyp+1,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
+ nterm_kcc(j,i)=2
+ nterm_kcc_Tb(j,i)=3
+ do k=1,nterm_kcc_Tb(j,i)
+ do l=1,nterm_kcc_Tb(j,i)
+ v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
+ +bnew1tor(k,2,i)*bnew2tor(l,2,j)
+ v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
+ +bnew1tor(k,2,i)*bnew2tor(l,1,j)
+ enddo
+ enddo
+ do k=1,nterm_kcc_Tb(j,i)
+ do l=1,nterm_kcc_Tb(j,i)
+#ifdef CORRCD
+ v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
+ -ccnewtor(k,2,i)*ddnewtor(l,2,j))
+ v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
+ +ccnewtor(k,1,i)*ddnewtor(l,2,j))
+#else
+ v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
+ -ccnewtor(k,2,i)*ddnewtor(l,2,j))
+ v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
+ +ccnewtor(k,1,i)*ddnewtor(l,2,j))
+#endif
+ enddo
+ enddo
+!c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
+ enddo
+ enddo
+#else
+ write (iout,*) "TOR_MODE>1 only with NEWCORR"
+ stop
+#endif
+ ENDIF ! TOR_MODE
+ if (tor_mode.gt.0 .and. lprint) then
+!c Print valence-torsional parameters
+ write (iout,'(a)') &
+ "Parameters of the valence-torsional potentials"
+ do i=-ntortyp+1,ntortyp-1
+ do j=-ntortyp+1,ntortyp-1
+ write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
+ write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
+ do k=1,nterm_kcc(j,i)
+ do l=1,nterm_kcc_Tb(j,i)
+ do ll=1,nterm_kcc_Tb(j,i)
+ write (iout,'(3i5,2f15.4)')&
+ k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ endif
+
#endif
!elwrite(iout,*) "parmread kontrol sc-bb"
! Read of Side-chain backbone correlation parameters
enddo
enddo
endif
-!
-! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
-! interaction energy of the Gly, Ala, and Pro prototypes.
-!
- read (ifourier,*) nloctyp
-!el 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(b1tilde(2,-nloctyp-1:nloctyp+1)) !(2,-maxtor:maxtor)
- allocate(cc(2,2,-nloctyp-1:nloctyp+1))
- allocate(dd(2,2,-nloctyp-1:nloctyp+1))
- 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)
- do i=1,2
- do ii=-nloctyp-1,nloctyp+1
- b1(i,ii)=0.0d0
- b2(i,ii)=0.0d0
- b1tilde(i,ii)=0.0d0
- do j=1,2
- cc(j,i,ii)=0.0d0
- dd(j,i,ii)=0.0d0
- ee(j,i,ii)=0.0d0
- ctilde(j,i,ii)=0.0d0
- dtilde(j,i,ii)=0.0d0
- enddo
- enddo
- enddo
-!--------------------------------
- allocate(b(13,0:nloctyp))
- do i=0,nloctyp-1
- read (ifourier,*)
- read (ifourier,*) (b(ii,i),ii=1,13)
- if (lprint) then
- write (iout,*) 'Type',i
- write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
- endif
- B1(1,i) = b(3,i)
- B1(2,i) = b(5,i)
- B1(1,-i) = b(3,i)
- B1(2,-i) = -b(5,i)
-! b1(1,i)=0.0d0
-! b1(2,i)=0.0d0
- B1tilde(1,i) = b(3,i)
- B1tilde(2,i) =-b(5,i)
- B1tilde(1,-i) =-b(3,i)
- B1tilde(2,-i) =b(5,i)
-! b1tilde(1,i)=0.0d0
-! b1tilde(2,i)=0.0d0
- B2(1,i) = b(2,i)
- B2(2,i) = b(4,i)
- B2(1,-i) =b(2,i)
- B2(2,-i) =-b(4,i)
-
-! b2(1,i)=0.0d0
-! b2(2,i)=0.0d0
- CC(1,1,i)= b(7,i)
- CC(2,2,i)=-b(7,i)
- CC(2,1,i)= b(9,i)
- CC(1,2,i)= b(9,i)
- CC(1,1,-i)= b(7,i)
- CC(2,2,-i)=-b(7,i)
- CC(2,1,-i)=-b(9,i)
- CC(1,2,-i)=-b(9,i)
-! CC(1,1,i)=0.0d0
-! CC(2,2,i)=0.0d0
-! CC(2,1,i)=0.0d0
-! CC(1,2,i)=0.0d0
- Ctilde(1,1,i)=b(7,i)
- Ctilde(1,2,i)=b(9,i)
- Ctilde(2,1,i)=-b(9,i)
- Ctilde(2,2,i)=b(7,i)
- Ctilde(1,1,-i)=b(7,i)
- Ctilde(1,2,-i)=-b(9,i)
- Ctilde(2,1,-i)=b(9,i)
- Ctilde(2,2,-i)=b(7,i)
-
-! Ctilde(1,1,i)=0.0d0
-! Ctilde(1,2,i)=0.0d0
-! Ctilde(2,1,i)=0.0d0
-! Ctilde(2,2,i)=0.0d0
- DD(1,1,i)= b(6,i)
- DD(2,2,i)=-b(6,i)
- DD(2,1,i)= b(8,i)
- DD(1,2,i)= b(8,i)
- DD(1,1,-i)= b(6,i)
- DD(2,2,-i)=-b(6,i)
- DD(2,1,-i)=-b(8,i)
- DD(1,2,-i)=-b(8,i)
-! DD(1,1,i)=0.0d0
-! DD(2,2,i)=0.0d0
-! DD(2,1,i)=0.0d0
-! DD(1,2,i)=0.0d0
- Dtilde(1,1,i)=b(6,i)
- Dtilde(1,2,i)=b(8,i)
- Dtilde(2,1,i)=-b(8,i)
- Dtilde(2,2,i)=b(6,i)
- Dtilde(1,1,-i)=b(6,i)
- Dtilde(1,2,-i)=-b(8,i)
- Dtilde(2,1,-i)=b(8,i)
- Dtilde(2,2,-i)=b(6,i)
-
-! Dtilde(1,1,i)=0.0d0
-! Dtilde(1,2,i)=0.0d0
-! Dtilde(2,1,i)=0.0d0
-! Dtilde(2,2,i)=0.0d0
- EE(1,1,i)= b(10,i)+b(11,i)
- EE(2,2,i)=-b(10,i)+b(11,i)
- EE(2,1,i)= b(12,i)-b(13,i)
- EE(1,2,i)= b(12,i)+b(13,i)
- EE(1,1,-i)= b(10,i)+b(11,i)
- EE(2,2,-i)=-b(10,i)+b(11,i)
- EE(2,1,-i)=-b(12,i)+b(13,i)
- EE(1,2,-i)=-b(12,i)-b(13,i)
-
-! ee(1,1,i)=1.0d0
-! ee(2,2,i)=1.0d0
-! ee(2,1,i)=0.0d0
-! ee(1,2,i)=0.0d0
-! ee(2,1,i)=ee(1,2,i)
-
- enddo
- if (lprint) then
- do i=1,nloctyp
- write (iout,*) 'Type',i
- write (iout,*) 'B1'
-! write (iout,'(f10.5)') B1(:,i)
- write(iout,*) B1(1,i),B1(2,i)
- write (iout,*) 'B2'
-! write (iout,'(f10.5)') B2(:,i)
- write(iout,*) B2(1,i),B2(2,i)
- write (iout,*) 'CC'
- do j=1,2
- write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
- enddo
- write(iout,*) 'DD'
- do j=1,2
- write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
- enddo
- write(iout,*) 'EE'
- do j=1,2
- write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
- enddo
- enddo
- endif
!
! Read electrostatic-interaction parameters
!
!---------------------- GB or BP potential -----------------------------
case (3:4)
! 30 do i=1,ntyp
+ if (scelemode.eq.0) then
do i=1,ntyp
read (isidep,*)(eps(i,j),j=i,ntyp)
enddo
write (iout,'(a3,6x,4f10.5)') (restyp(i,molnum(i)),sigma0(i),sigii(i),&
chip(i),alp(i),i=1,ntyp)
endif
+ else
+ allocate(icharge(ntyp1))
+! print *,ntyp,icharge(i)
+ icharge(:)=0
+ read (isidep,*) (icharge(i),i=1,ntyp)
+ print *,ntyp,icharge(i)
+! if(.not.allocated(eps)) allocate(eps(-ntyp
+ write (2,*) "icharge",(icharge(i),i=1,ntyp)
+ allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
+ allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
+ allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
+ allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
+ allocate(epsintab(ntyp,ntyp))
+ allocate(dtail(2,ntyp,ntyp))
+ print *,"control line 1"
+ allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
+ allocate(wqdip(2,ntyp,ntyp))
+ allocate(wstate(4,ntyp,ntyp))
+ allocate(dhead(2,2,ntyp,ntyp))
+ allocate(nstate(ntyp,ntyp))
+ allocate(debaykap(ntyp,ntyp))
+ print *,"control line 2"
+ if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
+ if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
+
+ do i=1,ntyp
+ do j=1,i
+! write (*,*) "Im in ALAB", i, " ", j
+ read(isidep,*) &
+ eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
+ (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
+ chis(i,j),chis(j,i), &
+ nstate(i,j),(wstate(k,i,j),k=1,4), &
+ dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
+ dtail(1,i,j),dtail(2,i,j), &
+ epshead(i,j),sig0head(i,j), &
+ rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
+ alphapol(i,j),alphapol(j,i), &
+ (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j)
+! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
+ END DO
+ END DO
+ DO i = 1, ntyp
+ DO j = i+1, ntyp
+ eps(i,j) = eps(j,i)
+ sigma(i,j) = sigma(j,i)
+ sigmap1(i,j)=sigmap1(j,i)
+ sigmap2(i,j)=sigmap2(j,i)
+ sigiso1(i,j)=sigiso1(j,i)
+ sigiso2(i,j)=sigiso2(j,i)
+! print *,"ATU",sigma(j,i),sigma(i,j),i,j
+ nstate(i,j) = nstate(j,i)
+ dtail(1,i,j) = dtail(1,j,i)
+ dtail(2,i,j) = dtail(2,j,i)
+ DO k = 1, 4
+ alphasur(k,i,j) = alphasur(k,j,i)
+ wstate(k,i,j) = wstate(k,j,i)
+ alphiso(k,i,j) = alphiso(k,j,i)
+ END DO
+
+ dhead(2,1,i,j) = dhead(1,1,j,i)
+ dhead(2,2,i,j) = dhead(1,2,j,i)
+ dhead(1,1,i,j) = dhead(2,1,j,i)
+ dhead(1,2,i,j) = dhead(2,2,j,i)
+
+ epshead(i,j) = epshead(j,i)
+ sig0head(i,j) = sig0head(j,i)
+
+ DO k = 1, 2
+ wqdip(k,i,j) = wqdip(k,j,i)
+ END DO
+
+ wquad(i,j) = wquad(j,i)
+ epsintab(i,j) = epsintab(j,i)
+ debaykap(i,j)=debaykap(j,i)
+! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
+ END DO
+ END DO
+ endif
+
! goto 50
!--------------------- GBV potential -----------------------------------
case (5)
!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)
+ if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1))
+ allocate(aa_aq(ntyp1,ntyp1),bb_aq(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)
+ if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
+ allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
+ allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
+ dcavtub(ntyp1))
+ allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
+ tubetranene(ntyp1))
do i=1,ntyp1
do j=1,ntyp1
aa_aq(i,j)=0.0D0
bb_aq(i,j)=0.0D0
aa_lip(i,j)=0.0d0
bb_lip(i,j)=0.0d0
+ if (scelemode.eq.0) then
chi(i,j)=0.0D0
sigma(i,j)=0.0D0
r0(i,j)=0.0D0
+ endif
enddo
enddo
!--------------------------------
epslip(i,j)=epslip(j,i)
enddo
enddo
+ if (scelemode.eq.0) then
do i=1,ntyp
do j=i,ntyp
sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
rs0(j,i)=rs0(i,j)
enddo
enddo
+ endif
if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
'Working parameters of the SC interactions:',&
' a ',' b ',' augm ',' sigma ',' r0 ',&
epsij=eps(i,j)
if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
rrij=sigma(i,j)
+ print *,"rrij",rrij
else
rrij=rr0(i)+rr0(j)
endif
bb_lip(i,j)=-sigeps*epsijlip*rrij
aa_lip(j,i)=aa_lip(i,j)
bb_lip(j,i)=bb_lip(i,j)
- if (ipot.gt.2) then
+ if ((ipot.gt.2).and. (scelemode.eq.0))then
sigt1sq=sigma0(i)**2
sigt2sq=sigma0(j)**2
sigii1=sigii(i)
endif
enddo
enddo
+ allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
+ allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
+ allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
+ allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
+ allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
+ allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
+ allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
+ allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
+
+! augm(:,:)=0.0D0
+! chip(:)=0.0D0
+! alp(:)=0.0D0
+! sigma0(:)=0.0D0
+! sigii(:)=0.0D0
+! rr0(:)=0.0D0
+
+ read (isidep_nucl,*) ipot_nucl
+! print *,"TU?!",ipot_nucl
+ if (ipot_nucl.eq.1) then
+ do i=1,ntyp_molec(2)
+ do j=i,ntyp_molec(2)
+ read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
+ elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
+ enddo
+ enddo
+ else
+ do i=1,ntyp_molec(2)
+ do j=i,ntyp_molec(2)
+ read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
+ chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
+ elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
+ enddo
+ enddo
+ endif
+! rpp(1,1)=2**(1.0/6.0)*5.16158
+ do i=1,ntyp_molec(2)
+ do j=i,ntyp_molec(2)
+ rrij=sigma_nucl(i,j)
+ r0_nucl(i,j)=rrij
+ r0_nucl(j,i)=rrij
+ rrij=rrij**expon
+ epsij=4*eps_nucl(i,j)
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ aa_nucl(i,j)=epsij*rrij*rrij
+ bb_nucl(i,j)=-sigeps*epsij*rrij
+ ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
+ ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
+ ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
+ ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
+ sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
+ 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
+ enddo
+ do j=1,i-1
+ aa_nucl(i,j)=aa_nucl(j,i)
+ bb_nucl(i,j)=bb_nucl(j,i)
+ ael3_nucl(i,j)=ael3_nucl(j,i)
+ ael6_nucl(i,j)=ael6_nucl(j,i)
+ ael63_nucl(i,j)=ael63_nucl(j,i)
+ ael32_nucl(i,j)=ael32_nucl(j,i)
+ elpp3_nucl(i,j)=elpp3_nucl(j,i)
+ elpp6_nucl(i,j)=elpp6_nucl(j,i)
+ elpp63_nucl(i,j)=elpp63_nucl(j,i)
+ elpp32_nucl(i,j)=elpp32_nucl(j,i)
+ eps_nucl(i,j)=eps_nucl(j,i)
+ sigma_nucl(i,j)=sigma_nucl(j,i)
+ sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
+ enddo
+ enddo
+
+ write(iout,*) "tube param"
+ read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
+ ccavtubpep,dcavtubpep,tubetranenepep
+ sigmapeptube=sigmapeptube**6
+ sigeps=dsign(1.0D0,epspeptube)
+ epspeptube=dabs(epspeptube)
+ pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
+ pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
+ write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
+ do i=1,ntyp
+ read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
+ ccavtub(i),dcavtub(i),tubetranene(i)
+ sigmasctube=sigmasctube**6
+ sigeps=dsign(1.0D0,epssctube)
+ epssctube=dabs(epssctube)
+ sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
+ sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
+ write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
+ enddo
+!-----------------READING SC BASE POTENTIALS-----------------------------
+ allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
+ allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
+ allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
+ allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
+ allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
+ allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
+ allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
+ allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
+ allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
+ allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
+ allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
+ allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
+ allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
+ allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
+ allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
+ allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
+
+
+ do i=1,ntyp_molec(1)
+ do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
+ write (*,*) "Im in ", i, " ", j
+ read(isidep_scbase,*) &
+ eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
+ chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
+ write(*,*) "eps",eps_scbase(i,j)
+ read(isidep_scbase,*) &
+ (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
+ chis_scbase(i,j,1),chis_scbase(i,j,2)
+ read(isidep_scbase,*) &
+ dhead_scbasei(i,j), &
+ dhead_scbasej(i,j), &
+ rborn_scbasei(i,j),rborn_scbasej(i,j)
+ read(isidep_scbase,*) &
+ (wdipdip_scbase(k,i,j),k=1,3), &
+ (wqdip_scbase(k,i,j),k=1,2)
+ read(isidep_scbase,*) &
+ alphapol_scbase(i,j), &
+ epsintab_scbase(i,j)
+ END DO
+ END DO
+ allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
+ allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
+
+ do i=1,ntyp_molec(1)
+ do j=1,ntyp_molec(2)-1
+ epsij=eps_scbase(i,j)
+ rrij=sigma_scbase(i,j)
+! r0(i,j)=rrij
+! r0(j,i)=rrij
+ rrij=rrij**expon
+! epsij=eps(i,j)
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ aa_scbase(i,j)=epsij*rrij*rrij
+ bb_scbase(i,j)=-sigeps*epsij*rrij
+ enddo
+ enddo
+!-----------------READING PEP BASE POTENTIALS-------------------
+ allocate(eps_pepbase(ntyp_molec(2)))
+ allocate(sigma_pepbase(ntyp_molec(2)))
+ allocate(chi_pepbase(ntyp_molec(2),2))
+ allocate(chipp_pepbase(ntyp_molec(2),2))
+ allocate(alphasur_pepbase(4,ntyp_molec(2)))
+ allocate(sigmap1_pepbase(ntyp_molec(2)))
+ allocate(sigmap2_pepbase(ntyp_molec(2)))
+ allocate(chis_pepbase(ntyp_molec(2),2))
+ allocate(wdipdip_pepbase(3,ntyp_molec(2)))
+
+
+ do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
+ write (*,*) "Im in ", i, " ", j
+ read(isidep_pepbase,*) &
+ eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
+ chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
+ write(*,*) "eps",eps_pepbase(j)
+ read(isidep_pepbase,*) &
+ (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
+ chis_pepbase(j,1),chis_pepbase(j,2)
+ read(isidep_pepbase,*) &
+ (wdipdip_pepbase(k,j),k=1,3)
+ END DO
+ allocate(aa_pepbase(ntyp_molec(2)))
+ allocate(bb_pepbase(ntyp_molec(2)))
+
+ do j=1,ntyp_molec(2)-1
+ epsij=eps_pepbase(j)
+ rrij=sigma_pepbase(j)
+! r0(i,j)=rrij
+! r0(j,i)=rrij
+ rrij=rrij**expon
+! epsij=eps(i,j)
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ aa_pepbase(j)=epsij*rrij*rrij
+ bb_pepbase(j)=-sigeps*epsij*rrij
+ enddo
+!--------------READING SC PHOSPHATE-------------------------------------
+!--------------READING SC PHOSPHATE-------------------------------------
+ allocate(eps_scpho(ntyp_molec(1)))
+ allocate(sigma_scpho(ntyp_molec(1)))
+ allocate(chi_scpho(ntyp_molec(1),2))
+ allocate(chipp_scpho(ntyp_molec(1),2))
+ allocate(alphasur_scpho(4,ntyp_molec(1)))
+ allocate(sigmap1_scpho(ntyp_molec(1)))
+ allocate(sigmap2_scpho(ntyp_molec(1)))
+ allocate(chis_scpho(ntyp_molec(1),2))
+ allocate(wqq_scpho(ntyp_molec(1)))
+ allocate(wqdip_scpho(2,ntyp_molec(1)))
+ allocate(alphapol_scpho(ntyp_molec(1)))
+ allocate(epsintab_scpho(ntyp_molec(1)))
+ allocate(dhead_scphoi(ntyp_molec(1)))
+ allocate(rborn_scphoi(ntyp_molec(1)))
+ allocate(rborn_scphoj(ntyp_molec(1)))
+ allocate(alphi_scpho(ntyp_molec(1)))
+
+
+! j=1
+ do j=1,ntyp_molec(1) ! without U then we will take T for U
+ write (*,*) "Im in scpho ", i, " ", j
+ read(isidep_scpho,*) &
+ eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
+ chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
+ write(*,*) "eps",eps_scpho(j)
+ read(isidep_scpho,*) &
+ (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
+ chis_scpho(j,1),chis_scpho(j,2)
+ read(isidep_scpho,*) &
+ (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
+ read(isidep_scpho,*) &
+ epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
+ alphi_scpho(j)
+
+ END DO
+ allocate(aa_scpho(ntyp_molec(1)))
+ allocate(bb_scpho(ntyp_molec(1)))
+ do j=1,ntyp_molec(1)
+ epsij=eps_scpho(j)
+ rrij=sigma_scpho(j)
+! r0(i,j)=rrij
+! r0(j,i)=rrij
+ rrij=rrij**expon
+! epsij=eps(i,j)
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ aa_scpho(j)=epsij*rrij*rrij
+ bb_scpho(j)=-sigeps*epsij*rrij
+ enddo
+
+
+ read(isidep_peppho,*) &
+ eps_peppho,sigma_peppho
+ read(isidep_peppho,*) &
+ (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
+ read(isidep_peppho,*) &
+ (wqdip_peppho(k),k=1,2)
+
+ epsij=eps_peppho
+ rrij=sigma_peppho
+! r0(i,j)=rrij
+! r0(j,i)=rrij
+ rrij=rrij**expon
+! epsij=eps(i,j)
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ aa_peppho=epsij*rrij*rrij
+ bb_peppho=-sigeps*epsij*rrij
+
allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
do i=1,ntyp
enddo
enddo
#endif
+ allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
+ read (itorp_nucl,*) ntortyp_nucl
+! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
+!el from energy module---------
+ allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
+ allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
+
+ allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
+ allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
+ allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
+ allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
+
+ allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
+ allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
+!el---------------------------
+ nterm_nucl(:,:)=0
+ nlor_nucl(:,:)=0
+!el--------------------
+ read (itorp_nucl,*) &
+ (itortyp_nucl(i),i=1,ntyp_molec(2))
+! print *,itortyp_nucl(:)
+!c write (iout,*) 'ntortyp',ntortyp
+ do i=1,ntortyp_nucl
+ do j=1,ntortyp_nucl
+ read (itorp_nucl,*) nterm_nucl(i,j),nlor_nucl(i,j)
+! print *,nterm_nucl(i,j),nlor_nucl(i,j)
+ v0ij=0.0d0
+ si=-1.0d0
+ do k=1,nterm_nucl(i,j)
+ read (itorp_nucl,*) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
+ v0ij=v0ij+si*v1_nucl(k,i,j)
+ si=-si
+ enddo
+ do k=1,nlor_nucl(i,j)
+ read (itorp_nucl,*) kk,vlor1_nucl(k,i,j),&
+ vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
+ v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
+ enddo
+ v0_nucl(i,j)=v0ij
+ enddo
+ enddo
+
!elwrite(iout,*) "parmread kontrol before oldscp"
!
enddo
endif
#endif
+ allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
+
+ do i=1,ntyp_molec(2)
+ read (iscpp_nucl,*) eps_scp_nucl(i),rscp_nucl(i)
+ enddo
+ do i=1,ntyp_molec(2)
+ aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
+ bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
+ enddo
+ r0pp=1.12246204830937298142*5.16158
+ epspp=4.95713/4
+ AEES=108.661
+ BEES=0.433246
+
!
! Define the constants of the disulfide bridge
!
!-----------------------------------------------------------------------------
subroutine read_general_data(*)
- use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions
- use energy_data, only:distchainmax
- use geometry_data, only:boxxsize,boxysize,boxzsize
+ use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,&
+ scelemode,TUBEmode,tor_mode
+
+ use energy_data, only:distchainmax,tubeR0,tubecenter
+ use geometry_data, only:boxxsize,boxysize,boxzsize,bordtubetop,&
+ bordtubebot,tubebufthick,buftubebot,buftubetop
! implicit none
! include "DIMENSIONS"
! include "DIMENSIONS.ZSCOPT"
call reada(controlcard,'BOXX',boxxsize,100.0d0)
call reada(controlcard,'BOXY',boxysize,100.0d0)
call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+ call readi(controlcard,"SCELEMODE",scelemode,0)
+ print *,"SCELE",scelemode
+ call readi(controlcard,'TORMODE',tor_mode,0)
+!C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
+ write(iout,*) "torsional and valence angle mode",tor_mode
+
+ call readi(controlcard,'TUBEMOD',tubemode,0)
+
+ if (TUBEmode.gt.0) then
+ call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
+ call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
+ call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
+ call reada(controlcard,"RTUBE",tubeR0,0.0d0)
+ call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
+ call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
+ call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
+ buftubebot=bordtubebot+tubebufthick
+ buftubetop=bordtubetop-tubebufthick
+ endif
ions=index(controlcard,"IONS").gt.0
call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
do k=1,3
cref(k,j+i,kkk)=cref(k,j,kkk)
enddo
+ write(iout,*) "J",j,"J+I",j+i
phi_ref(j+i)=phi_ref(j)
theta_ref(j+i)=theta_ref(j)
alph_ref(j+i)=alph_ref(j)
! include 'COMMON.HEADER'
! include 'COMMON.SBRIDGE'
character(len=50) :: tytul
- character(len=1),dimension(10) :: chainid=reshape((/'A','B','C',&
- 'D','E','F','G','H','I','J'/),shape(chainid))
+ character(len=1),dimension(24) :: chainid=reshape((/'A','B','C',&
+ 'D','E','F','G','H','I','J','K','L','M','N','O',&
+ 'P','Q','R','S','V','W','X','Y','Z'/),shape(chainid))
integer,dimension(nres) :: ica !(maxres)
real(kind=8) :: temp,efree,etot,entropy,rmsdev
integer :: ii,i,j,iti,ires,iatom,ichain,mnum
mnum=molnum(i)
iti=itype(i,mnum)
if (iti.eq.ntyp1) then
+ if (itype(i-1,molnum(i-1)).eq.ntyp1) then
ichain=ichain+1
ires=0
write (ipdb,'(a)') 'TER'
+ endif
else
ires=ires+1
iatom=iatom+1
projects / unres4.git / blobdiff