TMscore_subroutine.f
together.F
unres_csa.F
+ gnmr1.f
)
set(UNRES_CSA_SRC3 energy_p_new_barrier.F gradient_p.F )
& ns,nss,nfree,iss(maxss)
double precision dhpb,forcon
integer ihpb,jhpb,nhpb
- common /links/ dhpb(maxdim),forcon(maxdim),ihpb(maxdim),
- & jhpb(maxdim),nhpb
+ common /links/ dhpb(maxdim),dhpb1(maxdim),forcon(maxdim),
+ & ihpb(maxdim),jhpb(maxdim),ibecarb(maxdim),nhpb
double precision weidis
common /restraints/ weidis
integer link_start,link_end
together.o csa.o minim_jlee.o shift.o diff12.o bank.o newconf.o ran.o \
indexx.o prng_32.o contact.o gen_rand_conf.o \
sc_move.o test.o local_move.o rmsd.o fitsq.o elecont.o djacob.o \
- distfit.o banach.o TMscore_subroutine.o minim_mult.o refsys.o
+ distfit.o banach.o TMscore_subroutine.o minim_mult.o refsys.o gnmr1.o
no_option:
together.o csa.o minim_jlee.o shift.o diff12.o bank.o newconf.o ran.o \
indexx.o prng_32.o contact.o gen_rand_conf.o \
sc_move.o test.o local_move.o rmsd.o fitsq.o elecont.o djacob.o \
- distfit.o banach.o TMscore_subroutine.o minim_mult.o refsys.o
+ distfit.o banach.o TMscore_subroutine.o minim_mult.o refsys.o gnmr1.o
no_option:
together.o csa.o minim_jlee.o shift.o diff12.o bank.o newconf.o ran.o \
indexx.o prng_32.o contact.o gen_rand_conf.o \
sc_move.o test.o local_move.o rmsd.o fitsq.o elecont.o djacob.o \
- distfit.o banach.o TMscore_subroutine.o minim_mult.o refsys.o
+ distfit.o banach.o TMscore_subroutine.o minim_mult.o refsys.o gnmr1.o
no_option:
iii=ii
jjj=jj
endif
-cd write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj
+c write (iout,*) "i",i," ii",ii," iii",iii," jj",jj," jjj",jjj,
+c & dhpb(i),dhpb1(i),forcon(i)
C 24/11/03 AL: SS bridges handled separately because of introducing a specific
C distance and angle dependent SS bond potential.
- if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+cmc if (ii.gt.nres .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
+C 18/07/06 MC: Use the convention that the first nss pairs are SS bonds
+ if (i.le.nss) then
+ if (ii.gt.nres
+ & .and. itype(iii).eq.1 .and. itype(jjj).eq.1) then
call ssbond_ene(iii,jjj,eij)
ehpb=ehpb+2*eij
+ endif
cd write (iout,*) "eij",eij
+ else if (ii.gt.nres .and. jj.gt.nres) then
+c Restraints from contact prediction
+ dd=dist(ii,jj)
+ if (dhpb1(i).gt.0.0d0) then
+ ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c write (iout,*) "beta nmr",
+c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ else
+ dd=dist(ii,jj)
+ rdis=dd-dhpb(i)
+C Get the force constant corresponding to this distance.
+ waga=forcon(i)
+C Calculate the contribution to energy.
+ ehpb=ehpb+waga*rdis*rdis
+c write (iout,*) "beta reg",dd,waga*rdis*rdis
+C
+C Evaluate gradient.
+C
+ fac=waga*rdis/dd
+ endif
+ do j=1,3
+ ggg(j)=fac*(c(j,jj)-c(j,ii))
+ enddo
+ do j=1,3
+ ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+ ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+ enddo
+ do k=1,3
+ ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+ ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+ enddo
else
C Calculate the distance between the two points and its difference from the
C target distance.
dd=dist(ii,jj)
+ if (dhpb1(i).gt.0.0d0) then
+ ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd
+c write (iout,*) "alph nmr",
+c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ else
rdis=dd-dhpb(i)
C Get the force constant corresponding to this distance.
waga=forcon(i)
C Calculate the contribution to energy.
ehpb=ehpb+waga*rdis*rdis
+c write (iout,*) "alpha reg",dd,waga*rdis*rdis
C
C Evaluate gradient.
C
fac=waga*rdis/dd
+ endif
cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
cd & ' waga=',waga,' fac=',fac
do j=1,3
--- /dev/null
+ double precision function gnmr1(y,ymin,ymax)
+ implicit none
+ double precision y,ymin,ymax
+ double precision wykl /4.0d0/
+ if (y.lt.ymin) then
+ gnmr1=(ymin-y)**wykl/wykl
+ else if (y.gt.ymax) then
+ gnmr1=(y-ymax)**wykl/wykl
+ else
+ gnmr1=0.0d0
+ endif
+ return
+ end
+c------------------------------------------------------------------------------
+ double precision function gnmr1prim(y,ymin,ymax)
+ implicit none
+ double precision y,ymin,ymax
+ double precision wykl /4.0d0/
+ if (y.lt.ymin) then
+ gnmr1prim=-(ymin-y)**(wykl-1)
+ else if (y.gt.ymax) then
+ gnmr1prim=(y-ymax)**(wykl-1)
+ else
+ gnmr1prim=0.0d0
+ endif
+ return
+ end
+c------------------------------------------------------------------------------
+ double precision function harmonic(y,ymax)
+ implicit none
+ double precision y,ymax
+ double precision wykl /2.0d0/
+ harmonic=(y-ymax)**wykl
+ return
+ end
+c-------------------------------------------------------------------------------
+ double precision function harmonicprim(y,ymax)
+ double precision y,ymin,ymax
+ double precision wykl /2.0d0/
+ harmonicprim=(y-ymax)*wykl
+ return
+ end
+c---------------------------------------------------------------------------------
enddo
call contact(.true.,ncont_ref,icont_ref,co)
endif
-c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
- call flush(iout)
- if (constr_dist.gt.0) call read_dist_constr
-c write (iout,*) "After read_dist_constr nhpb",nhpb
- call hpb_partition
if(me.eq.king.or..not.out1file)
& write (iout,*) 'Contact order:',co
if (pdbref) then
enddo
endif
endif
+c write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
+ if (constr_dist.gt.0) then
+ call read_dist_constr
+ endif
+ if (nhpb.gt.0) call hpb_partition
+c write (iout,*) "After read_dist_constr nhpb",nhpb
+c call flush(iout)
if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
& .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
& modecalc.ne.10) then
integer ifrag_(2,100),ipair_(2,100)
double precision wfrag_(100),wpair_(100)
character*500 controlcard
- write (iout,*) "Calling read_dist_constr"
- write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
- call flush(iout)
+c write (iout,*) "Calling read_dist_constr"
+c write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
+c call flush(iout)
call card_concat(controlcard)
call readi(controlcard,"NFRAG",nfrag_,0)
call readi(controlcard,"NPAIR",npair_,0)
c do i=1,npair_
c write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
c enddo
+ if (.not.refstr .and. nfrag.gt.0) then
+ write (iout,*)
+ & "ERROR: no reference structure to compute distance restraints"
+ write (iout,*)
+ & "Restraints must be specified explicitly (NDIST=number)"
+ stop
+ endif
+ if (nfrag.lt.2 .and. npair.gt.0) then
+ write (iout,*) "ERROR: Less than 2 fragments specified",
+ & " but distance restraints between pairs requested"
+ stop
+ endif
call flush(iout)
do i=1,nfrag_
if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
if (wfrag_(i).gt.0.0d0) then
do j=ifrag_(1,i),ifrag_(2,i)-1
do k=j+1,ifrag_(2,i)
- write (iout,*) "j",j," k",k
+c write (iout,*) "j",j," k",k
ddjk=dist(j,k)
if (constr_dist.eq.1) then
nhpb=nhpb+1
endif
enddo
do i=1,ndist_
- read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
+ read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
+ & ibecarb(i),forcon(nhpb+1)
if (forcon(nhpb+1).gt.0.0d0) then
nhpb=nhpb+1
- dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+ if (ibecarb(i).gt.0) then
+ ihpb(i)=ihpb(i)+nres
+ jhpb(i)=jhpb(i)+nres
+ endif
+ if (dhpb(nhpb).eq.0.0d0)
+ & dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
+ endif
+ enddo
#ifdef MPI
- if (.not.out1file .or. me.eq.king)
- & write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
- & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
-#else
- write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
- & nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
+ if (.not.out1file .or. me.eq.king) then
#endif
- endif
+ do i=1,nhpb
+ write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
+ & i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
enddo
call flush(iout)
+#ifdef MPI
+ endif
+#endif
return
end
c-------------------------------------------------------------------------------
projects / unres.git / commitdiff