endif
enddo
#else
+ mp(:)=0.0d0
+ ip(:)=0.0d0
+ msc(:,:)=0.0d0
+ isc(:,:)=0.0d0
+
allocate(msc(ntyp+1,5)) !(ntyp+1)
allocate(isc(ntyp+1,5)) !(ntyp+1)
allocate(restok(ntyp+1,5)) !(ntyp+1)
!---------------------- GB or BP potential -----------------------------
case(3:4)
! 30 do i=1,ntyp
+! print *,"I AM in SCELE",scelemode
+ if (scelemode.eq.0) then
do i=1,ntyp
read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
enddo
write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
chip(i),alp(i),i=1,ntyp)
endif
+ else
+! print *,ntyp,"NTYP"
+ 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
+!c 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))
+ 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))
+ 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)
+! 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)
+! 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)
! Calculate the "working" parameters of SC interactions.
!el from module energy - COMMON.INTERACT-------
- allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1),chi(ntyp1,ntyp1)) !(ntyp,ntyp)
+ allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
+ if (.not.allocated(chi)) allocate(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)
+ 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),&
bb_aq(:,:)=0.0D0
aa_lip(:,:)=0.0D0
bb_lip(:,:)=0.0D0
+ if (scelemode.eq.0) then
chi(:,:)=0.0D0
sigma(:,:)=0.0D0
r0(:,:)=0.0D0
+ endif
acavtub(:)=0.0d0
bcavtub(:)=0.0d0
ccavtub(:)=0.0d0
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 *,"SIGMA ZLA?",sigma(i,j)
else
rrij=rr0(i)+rr0(j)
endif
sigeps=dsign(1.0D0,epsij)
epsij=dabs(epsij)
aa_aq(i,j)=epsij*rrij*rrij
+ print *,"ADASKO",epsij,rrij,expon
bb_aq(i,j)=-sigeps*epsij*rrij
aa_aq(j,i)=aa_aq(i,j)
bb_aq(j,i)=bb_aq(i,j)
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
+ if ((ipot.gt.2).and. (scelemode.eq.0)) then
sigt1sq=sigma0(i)**2
sigt2sq=sigma0(j)**2
sigii1=sigii(i)
! system
nres=ires
if (((ires_old.ne.ires).and.(molecule.ne.5)) &
- .and.(.not.unres_pdb)) &
+ ) &
nres_molec(molecule)=nres_molec(molecule)-2
print *,'I have',nres, nres_molec(:)
istype(i)=istype_temp(i)
enddo
enddo
- if (itype(1,1).eq.ntyp1) then
- nsup=nsup-1
- nstart_sup=2
- if (unres_pdb) then
+! if (itype(1,1).eq.ntyp1) then
+! nsup=nsup-1
+! nstart_sup=2
+! if (unres_pdb) then
! 2/15/2013 by Adam: corrected insertion of the first dummy residue
- call refsys(2,3,4,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,1)=c(j,2)-1.9d0*e2(j)
- enddo
- else
- do j=1,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
- endif
- endif
+! call refsys(2,3,4,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,1)=c(j,2)-1.9d0*e2(j)
+! enddo
+! else
+! do j=1,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
+! endif
+! endif
if (lprn) then
write (iout,'(/a)') &
call reada(controlcard,'BOXY',boxysize,100.0d0)
call reada(controlcard,'BOXZ',boxzsize,100.0d0)
call readi(controlcard,'TUBEMOD',tubemode,0)
+ print *,"SCELE",scelemode
+ call readi(controlcard,"SCELEMODE",scelemode,0)
+ print *,"SCELE",scelemode
+
+! elemode = 0 is orignal UNRES electrostatics
+! elemode = 1 is "Momo" potentials in progress
+! elemode = 2 is in development EVALD
write (iout,*) TUBEmode,"TUBEMODE"
if (TUBEmode.gt.0) then
call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
projects / unres4.git / blobdiff