call mygetenv('SCPPAR_NUCL',scpname_nucl)
open (iscpp_nucl,file=scpname_nucl,status='old')
+ call mygetenv('IONPAR_NUCL',ionnuclname)
+ open (iionnucl,file=ionnuclname,status='old')
#ifndef OLDSCP
!
use energy_data
use geometry_data, only:nres,deg2rad,c,dc,nres_molec
- use control_data, only:iscode
+ use control_data, only:iscode,dyn_ss
use io_base, only:rescode
- use control, only:setup_var,init_int_table
+ use control, only:setup_var,init_int_table,hpb_partition
use geometry, only:alloc_geo_arrays
use energy, only:alloc_ener_arrays
! implicit real*8 (a-h,o-z)
dhpb(i),ebr,forcon(i)
enddo
endif
+ if (ns.gt.0.and.dyn_ss) then
+ do i=nss+1,nhpb
+ ihpb(i-nss)=ihpb(i)
+ jhpb(i-nss)=jhpb(i)
+ forcon(i-nss)=forcon(i)
+ dhpb(i-nss)=dhpb(i)
+ enddo
+ nhpb=nhpb-nss
+ nss=0
+ call hpb_partition
+ do i=1,ns
+ dyn_ss_mask(iss(i))=.true.
+ enddo
+ endif
write (iout,'(a)')
return
end subroutine molread
real(kind=8) :: dwa16,akl,si,rri,epsij,rrij,sigeps,sigt1sq,epsijlip,&
sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm, &
epspeptube,epssctube,sigmapeptube, &
- sigmasctube
+ sigmasctube,ssscale
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,jj
integer :: nlobi,iblock,maxinter,iscprol,ncatprotparm
character*3 string
-
+
!
! Body
!
vbl=3.8D0
vblinv=1.0D0/vbl
vblinv2=vblinv*vblinv
+ itime_mat=0
#ifndef CLUSTER
call card_concat(controlcard,.true.)
wname(4)="WCORRH"
!el
+ call reada(controlcard,"D0CM",d0cm,3.78d0)
+ call reada(controlcard,"AKCM",akcm,15.1d0)
+ call reada(controlcard,"AKTH",akth,11.0d0)
+ call reada(controlcard,"AKCT",akct,12.0d0)
+ call reada(controlcard,"V1SS",v1ss,-1.08d0)
+ call reada(controlcard,"V2SS",v2ss,7.61d0)
+ call reada(controlcard,"V3SS",v3ss,13.7d0)
+ call reada(controlcard,"EBR",ebr,-5.50D0)
+ call reada(controlcard,"ATRISS",atriss,0.301D0)
+ call reada(controlcard,"BTRISS",btriss,0.021D0)
+ call reada(controlcard,"CTRISS",ctriss,1.001D0)
+ call reada(controlcard,"DTRISS",dtriss,1.001D0)
+ call reada(controlcard,"SSSCALE",ssscale,1.0D0)
+ dyn_ss=(index(controlcard,'DYN_SS').gt.0)
+
allocate(ww(max_eneW))
do i=1,n_eneW
key = wname(i)(:ilen(wname(i)))
wpepbase=ww(47)
wscpho=ww(48)
wpeppho=ww(49)
+ wcatnucl=ww(50)
! print *,"KURWA",ww(48)
! "WSCBASE ","WPEPBASE ","WSCPHO ","WPEPPHO "
! "WVDWPP ","WELPP ","WVDWPSB ","WELPSB ","WVDWSB ",&
weights(47)= wpepbase
weights(48) =wscpho
weights(49) =wpeppho
+ weights(50) =wcatnucl
+
! el--------
call card_concat(controlcard,.false.)
write (iout,*) "Parameter set:",iparm
write (iout,*) "Energy-term weights:"
do i=1,n_eneW
- write (iout,'(a16,f10.5)') wname(i),ww(i)
+ write (iout,'(i3,a16,f10.5)') i,wname(i),ww(i)
enddo
write (iout,*) "Sidechain potential file : ",&
sidename_t(:ilen(sidename_t))
endif
if (.not. allocated(msc)) allocate(msc(ntyp1,5))
if (.not. allocated(restok)) allocate(restok(ntyp1,5))
+ if (oldion.eq.1) then
do i=1,ntyp_molec(5)
read(iion,*) msc(i,5),restok(i,5)
read (iion,*) (catprm(i,k),i=1,ncatprotparm)
enddo
print *, catprm
+ endif
+ allocate(catnuclprm(14,ntyp_molec(2),ntyp_molec(5)))
+ do i=1,ntyp_molec(5)
+ do j=1,ntyp_molec(2)
+ write(iout,*) i,j
+ read(iionnucl,*) (catnuclprm(k,j,i),k=1,14)
+ enddo
+ enddo
+ write(*,'(3(5x,a6)11(7x,a6))') "w1 ","w2 ","epslj ","pis1 ", &
+ "sigma0","epsi0 ","chi1 ","chip1 ","sig ","b1 ","b2 ", &
+ "b3 ","b4 ","chis1 "
+ do i=1,ntyp_molec(5)
+ do j=1,ntyp_molec(2)
+ write(*,'(3(f10.3,x),11(f12.6,x),a3,2a)') (catnuclprm(k,j,i),k=1,14), &
+ restyp(i,5),"-",restyp(j,2)
+ enddo
+ enddo
+
read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
do i=1,ntyp_molec(2)
nbondterm_nucl(i)=1
bad(i,j)=0.0D0
enddo
enddo
+! Ions by Aga
+
+ allocate(alphapolcat(ntyp,ntyp),epsheadcat(ntyp,ntyp),sig0headcat(ntyp,ntyp))
+ allocate(alphapolcat2(ntyp,ntyp))
+ allocate(sigiso1cat(ntyp,ntyp),rborn1cat(ntyp,ntyp),rborn2cat(ntyp,ntyp),sigmap1cat(ntyp,ntyp))
+ allocate(sigmap2cat(ntyp,ntyp),sigiso2cat(ntyp,ntyp))
+ allocate(chis1cat(ntyp,ntyp),chis2cat(ntyp,ntyp),wquadcat(ntyp,ntyp),chipp1cat(ntyp,ntyp),chipp2cat(ntyp,ntyp))
+ allocate(epsintabcat(ntyp,ntyp))
+ allocate(dtailcat(2,ntyp,ntyp))
+ allocate(alphasurcat(4,ntyp,ntyp),alphisocat(4,ntyp,ntyp))
+ allocate(wqdipcat(2,ntyp,ntyp))
+ allocate(wstatecat(4,ntyp,ntyp))
+ allocate(dheadcat(2,2,ntyp,ntyp))
+ allocate(nstatecat(ntyp,ntyp))
+ allocate(debaykapcat(ntyp,ntyp))
+
+ if (.not.allocated(epscat)) allocate (epscat(0:ntyp1,0:ntyp1))
+ if (.not.allocated(sigmacat)) allocate(sigmacat(0:ntyp1,0:ntyp1))
+! if (.not.allocated(chicat)) allocate(chicat(ntyp1,ntyp1)) !(ntyp,ntyp)
+ if (.not.allocated(chi1cat)) allocate(chi1cat(ntyp1,ntyp1)) !(ntyp,ntyp)
+ if (.not.allocated(chi2cat)) allocate(chi2cat(ntyp1,ntyp1)) !(ntyp,ntyp)
+
+
+ allocate (ichargecat(ntyp_molec(5)))
+! i to SC, j to jon, isideocat - nazwa pliku z ktorego czytam parametry
+ if (oldion.eq.0) then
+ if (.not.allocated(icharge)) then ! this mean you are oprating in old sc-sc mode
+ allocate(icharge(1:ntyp1))
+ read(iion,*) (icharge(i),i=1,ntyp)
+ else
+ read(iion,*) ijunk
+ endif
+
+ do i=1,ntyp_molec(5)
+ read(iion,*) msc(i,5),restok(i,5),ichargecat(i)
+ print *,msc(i,5),restok(i,5)
+ enddo
+ ip(5)=0.2
+!DIR$ NOUNROLL
+ do j=1,ntyp_molec(5)
+ do i=1,ntyp
+! do j=1,ntyp_molec(5)
+! write (*,*) "Im in ALAB", i, " ", j
+ read(iion,*) &
+ epscat(i,j),sigmacat(i,j), &
+! chicat(i,j),chicat(j,i),chippcat(i,j),chippcat(j,i), &
+ chi1cat(i,j),chi2cat(i,j),chipp1cat(i,j),chipp2cat(i,j), &
+
+ (alphasurcat(k,i,j),k=1,4),sigmap1cat(i,j),sigmap2cat(i,j),&
+! chiscat(i,j),chiscat(j,i), &
+ chis1cat(i,j),chis2cat(i,j), &
+
+ nstatecat(i,j),(wstatecat(k,i,j),k=1,4), & !5 w tej lini - 1 integer pierwszy
+ dheadcat(1,1,i,j),dheadcat(1,2,i,j),dheadcat(2,1,i,j),dheadcat(2,2,i,j),&
+ dtailcat(1,i,j),dtailcat(2,i,j), &
+ epsheadcat(i,j),sig0headcat(i,j), &
+!wdipcat = w1 , w2
+! rborncat(i,j),rborncat(j,i),&
+ rborn1cat(i,j),rborn2cat(i,j),&
+ (wqdipcat(k,i,j),k=1,2), &
+ alphapolcat(i,j),alphapolcat2(j,i), &
+ (alphisocat(k,i,j),k=1,4),sigiso1cat(i,j),sigiso2cat(i,j),epsintabcat(i,j),debaykapcat(i,j)
+! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
+ END DO
+ END DO
+ allocate(aa_aq_cat(-ntyp:ntyp,ntyp),bb_aq_cat(-ntyp:ntyp,ntyp))
+ do i=1,ntyp
+ do j=1,ntyp_molec(5)
+ epsij=epscat(i,j)
+ rrij=sigmacat(i,j)
+ rrij=rrij**expon
+ sigeps=dsign(1.0D0,epsij)
+ epsij=dabs(epsij)
+ aa_aq_cat(i,j)=epsij*rrij*rrij
+ bb_aq_cat(i,j)=-sigeps*epsij*rrij
+ enddo
+ enddo
+ do i=1,ntyp
+ do j=1,ntyp_molec(5)
+ if (i.eq.10) then
+ write (iout,*) 'i= ', i, ' j= ', j
+ write (iout,*) 'epsi0= ', epscat(i,j)
+ write (iout,*) 'sigma0= ', sigmacat(i,j)
+ write (iout,*) 'chi1= ', chi1cat(i,j)
+ write (iout,*) 'chi1= ', chi2cat(i,j)
+ write (iout,*) 'chip1= ', chipp1cat(1,j)
+ write (iout,*) 'chip2= ', chipp2cat(1,j)
+ write (iout,*) 'alphasur1= ', alphasurcat(1,1,j)
+ write (iout,*) 'alphasur2= ', alphasurcat(2,1,j)
+ write (iout,*) 'alphasur3= ', alphasurcat(3,1,j)
+ write (iout,*) 'alphasur4= ', alphasurcat(4,1,j)
+ write (iout,*) 'sig1= ', sigmap1cat(1,j)
+ write (iout,*) 'sig2= ', sigmap2cat(1,j)
+ write (iout,*) 'chis1= ', chis1cat(1,j)
+ write (iout,*) 'chis1= ', chis2cat(1,j)
+ write (iout,*) 'nstatecat(i,j)= ', nstatecat(1,j)
+ write (iout,*) 'wstatecat(k,i,j),k=1= ',wstatecat(1,1,j)
+ write (iout,*) 'dhead= ', dheadcat(1,1,1,j)
+ write (iout,*) 'dhead2= ', dheadcat(1,2,1,j)
+ write (iout,*) 'a1= ', rborn1cat(i,j)
+ write (iout,*) 'a2= ', rborn2cat(i,j)
+ write (iout,*) 'epsin= ', epsintabcat(1,j), epsintabcat(j,1)
+ write (iout,*) 'alphapol1= ', alphapolcat(1,j)
+ write (iout,*) 'alphapol2= ', alphapolcat(j,1)
+ write (iout,*) 'w1= ', wqdipcat(1,i,j)
+ write (iout,*) 'w2= ', wqdipcat(2,i,j)
+ write (iout,*) 'debaykapcat(i,j)= ', debaykapcat(1,j)
+ endif
+
+ If ((i.eq.1).and.(j.eq.27)) then
+ write (iout,*) 'SEP'
+ Write (iout,*) 'w1= ', wqdipcat(1,1,27)
+ Write (iout,*) 'w2= ', wqdipcat(2,1,27)
+ endif
+
+ enddo
+ enddo
+
+ endif
+
#ifdef CLUSTER
!
! Define the SC-p interaction constants
!
! Define the constants of the disulfide bridge
!
- ebr=-5.50D0
+! ebr=-5.50D0
!
! Old arbitrary potential - commented out.
!
! energy surface of diethyl disulfide.
! A. Liwo and U. Kozlowska, 11/24/03
!
- D0CM = 3.78d0
- AKCM = 15.1d0
- AKTH = 11.0d0
- AKCT = 12.0d0
- V1SS =-1.08d0
- V2SS = 7.61d0
- V3SS = 13.7d0
+! D0CM = 3.78d0
+! AKCM = 15.1d0
+! AKTH = 11.0d0
+! AKCT = 12.0d0
+! V1SS =-1.08d0
+! V2SS = 7.61d0
+! V3SS = 13.7d0
+#ifndef CLUSTER
+ if (dyn_ss) then
+ ss_depth=(ebr/wsc-0.25*eps(1,1))*ssscale
+ Ht=(Ht/wsc-0.25*eps(1,1))*ssscale
+ akcm=akcm/wsc*ssscale
+ akth=akth/wsc*ssscale
+ akct=akct/wsc*ssscale
+ v1ss=v1ss/wsc*ssscale
+ v2ss=v2ss/wsc*ssscale
+ v3ss=v3ss/wsc*ssscale
+ else
+ ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
+ endif
+#endif
if (lprint) then
write (iout,'(/a)') "Disulfide bridge parameters:"
write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
subroutine read_general_data(*)
use control_data, only:indpdb,symetr,r_cut_ele,rlamb_ele,ions,&
- scelemode,TUBEmode,tor_mode
+ scelemode,TUBEmode,tor_mode,energy_dec
- use energy_data, only:distchainmax,tubeR0,tubecenter
+ use energy_data, only:distchainmax,tubeR0,tubecenter,dyn_ss
use geometry_data, only:boxxsize,boxysize,boxzsize,bordtubetop,&
bordtubebot,tubebufthick,buftubebot,buftubetop
! implicit none
call reada(controlcard,'BOXX',boxxsize,100.0d0)
call reada(controlcard,'BOXY',boxysize,100.0d0)
call reada(controlcard,'BOXZ',boxzsize,100.0d0)
+ energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
call readi(controlcard,"SCELEMODE",scelemode,0)
+ call readi(controlcard,"OLDION",oldion,0)
+ dyn_ss=(index(controlcard,'DYN_SS').gt.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
buftubetop=bordtubetop-tubebufthick
endif
ions=index(controlcard,"IONS").gt.0
- call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
+ call reada(controlcard,"R_CUT_ELE",r_cut_ele,25.0d0)
call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
write(iout,*) "R_CUT_ELE=",r_cut_ele
check_conf=index(controlcard,"NO_CHECK_CONF").eq.0
!--------------------------------------------------------------------------------
subroutine pdboutW(ii,temp,efree,etot,entropy,rmsdev)
- use geometry_data, only:nres,c
+ use geometry_data, only:nres,c,boxxsize,boxysize,boxzsize
use energy_data, only:nss,nnt,nct,ihpb,jhpb,itype,molnum
+ use energy, only:boxshift
! implicit real*8 (a-h,o-z)
! include 'DIMENSIONS'
! include 'DIMENSIONS.ZSCOPT'
'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
+ real(kind=8) :: temp,efree,etot,entropy,rmsdev,xj,yj,zj
integer :: ii,i,j,iti,ires,iatom,ichain,mnum
write(ipdb,'("REMARK CONF",i8," TEMPERATURE",f7.1," RMS",0pf7.2)')&
ii,temp,rmsdev
ires=ires+1
iatom=iatom+1
ica(i)=iatom
+ if (mnum.ne.5) then
write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
ires,(c(j,i),j=1,3)
- if (iti.ne.10) then
+ else
+ xj=boxshift(c(1,i)-c(1,2),boxxsize)
+ yj=boxshift(c(2,i)-c(2,2),boxysize)
+ zj=boxshift(c(3,i)-c(3,2),boxzsize)
+ write (ipdb,10) iatom,restyp(iti,mnum),chainid(ichain),&
+ ires,c(1,2)+xj,c(2,2)+yj,c(3,2)+zj
+ endif
+ if ((iti.ne.10).and.(mnum.ne.5)) then
iatom=iatom+1
write (ipdb,20) iatom,restyp(iti,mnum),chainid(ichain),&
ires,(c(j,nres+i),j=1,3)
projects / unres4.git / blobdiff