X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-NEWSC%2Fqwolynes.f;fp=source%2Fwham%2Fsrc-NEWSC%2Fqwolynes.f;h=97b5efbb535963c11c4a875300ad4206e8cb754a;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/wham/src-NEWSC/qwolynes.f b/source/wham/src-NEWSC/qwolynes.f new file mode 100755 index 0000000..97b5efb --- /dev/null +++ b/source/wham/src-NEWSC/qwolynes.f @@ -0,0 +1,186 @@ + double precision function qwolynes(ilevel,jfrag) + implicit none + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'DIMENSIONS.COMPAR' + include 'COMMON.IOUNITS' + include 'COMMON.COMPAR' + include 'COMMON.CHAIN' + include 'COMMON.INTERACT' + integer ilevel,jfrag + integer i,j,jl,k,l,il,kl,nl,np,ip,kp + integer nsep /3/ + double precision dist + double precision qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM + logical lprn /.false./ + double precision sigm,x + sigm(x)=0.25d0*x +c write (iout,*) "QWolyes: " jfrag",jfrag, +c & " ilevel",ilevel + qq = 0.0d0 + if (ilevel.eq.0) then + if (lprn) write (iout,*) "Q computed for whole molecule" + nl=0 + do il=nnt+nsep,nct + do jl=nnt,il-nsep + dij=0.0d0 + dijCM=0.0d0 + d0ij=0.0d0 + d0ijCM=0.0d0 + qqij=0.0d0 + qqijCM=0.0d0 + nl=nl+1 + d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+ + & (cref(2,jl)-cref(2,il))**2+ + & (cref(3,jl)-cref(3,il))**2) + dij=dist(il,jl) + qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2) + if (itype(il).ne.10 .or. itype(jl).ne.10) then + nl=nl+1 + d0ijCM=dsqrt( + & (cref(1,jl+nres)-cref(1,il+nres))**2+ + & (cref(2,jl+nres)-cref(2,il+nres))**2+ + & (cref(3,jl+nres)-cref(3,il+nres))**2) + dijCM=dist(il+nres,jl+nres) + qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2) + endif + qq = qq+qqij+qqijCM + if (lprn) then + write (iout,*) "il",il," jl",jl, + & " itype",itype(il),itype(jl) + write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM, + & " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM + endif + enddo + enddo + qq = qq/nl + if (lprn) write (iout,*) "nl",nl," qq",qq + else if (ilevel.eq.1) then + if (lprn) write (iout,*) "Level",ilevel," fragment",jfrag + nl=0 +c write (iout,*) "nlist_frag",nlist_frag(jfrag) + do i=2,nlist_frag(jfrag) + do j=1,i-1 + il=list_frag(i,jfrag) + jl=list_frag(j,jfrag) + if (iabs(il-jl).gt.nsep) then + dij=0.0d0 + dijCM=0.0d0 + d0ij=0.0d0 + d0ijCM=0.0d0 + qqij=0.0d0 + qqijCM=0.0d0 + nl=nl+1 + d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+ + & (cref(2,jl)-cref(2,il))**2+ + & (cref(3,jl)-cref(3,il))**2) + dij=dist(il,jl) + qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2) + if (itype(il).ne.10 .or. itype(jl).ne.10) then + nl=nl+1 + d0ijCM=dsqrt( + & (cref(1,jl+nres)-cref(1,il+nres))**2+ + & (cref(2,jl+nres)-cref(2,il+nres))**2+ + & (cref(3,jl+nres)-cref(3,il+nres))**2) + dijCM=dist(il+nres,jl+nres) + qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2) + endif + qq = qq+qqij+qqijCM + if (lprn) then + write (iout,*) "i",i," j",j," il",il," jl",jl, + & " itype",itype(il),itype(jl) + write (iout,*)"d0ij",d0ij," dij",dij," d0ijCM",d0ijCM, + & " dijCM",dijCM," qqij",qqij," qqijCM",qqijCM + endif + endif + enddo + enddo + qq = qq/nl + if (lprn) write (iout,*) "nl",nl," qq",qq + else if (ilevel.eq.2) then + np=npiece(jfrag,ilevel) + nl=0 + do i=2,np + ip=ipiece(i,jfrag,ilevel) + do j=1,nlist_frag(ip) + il=list_frag(j,ip) + do k=1,i-1 + kp=ipiece(k,jfrag,ilevel) + do l=1,nlist_frag(kp) + kl=list_frag(l,kp) + if (iabs(kl-il).gt.nsep) then + nl=nl+1 + dij=0.0d0 + dijCM=0.0d0 + d0ij=0.0d0 + d0ijCM=0.0d0 + qqij=0.0d0 + qqijCM=0.0d0 + d0ij=dsqrt((cref(1,kl)-cref(1,il))**2+ + & (cref(2,kl)-cref(2,il))**2+ + & (cref(3,kl)-cref(3,il))**2) + dij=dist(il,kl) + qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2) + if (itype(il).ne.10 .or. itype(kl).ne.10) then + nl=nl+1 + d0ijCM=dsqrt( + & (cref(1,kl+nres)-cref(1,il+nres))**2+ + & (cref(2,kl+nres)-cref(2,il+nres))**2+ + & (cref(3,kl+nres)-cref(3,il+nres))**2) + dijCM=dist(il+nres,kl+nres) + qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/ + & (sigm(d0ijCM)))**2) + endif + qq = qq+qqij+qqijCM + if (lprn) then + write (iout,*) "i",i," j",j," k",k," l",l," il",il, + & " kl",kl," itype",itype(il),itype(kl) + write (iout,*) " d0ij",d0ij," dij",dij," d0ijCM", + & d0ijCM," dijCM",dijCM," qqij",qqij," qqijCM",qqijCM + endif + endif + enddo ! l + enddo ! k + enddo ! j + enddo ! i + qq = qq/nl + if (lprn) write (iout,*) "nl",nl," qq",qq + else + write (iout,*)"Error: Q can be computed only for level 1 and 2." + endif + qwolynes=1.0d0-qq + return + end +c------------------------------------------------------------------------------- + subroutine fragment_list + implicit none + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'DIMENSIONS.COMPAR' + include 'COMMON.IOUNITS' + include 'COMMON.COMPAR' + logical lprn /.true./ + integer i,ilevel,j,k,jfrag + do jfrag=1,nfrag(1) + nlist_frag(jfrag)=0 + do i=1,npiece(jfrag,1) + if (lprn) write (iout,*) "jfrag=",jfrag, + & "i=",i," fragment",ifrag(1,i,jfrag), + & ifrag(2,i,jfrag) + do j=ifrag(1,i,jfrag),ifrag(2,i,jfrag) + do k=1,nlist_frag(jfrag) + if (list_frag(k,jfrag).eq.j) goto 10 + enddo + nlist_frag(jfrag)=nlist_frag(jfrag)+1 + list_frag(nlist_frag(jfrag),jfrag)=j + enddo + 10 continue + enddo + enddo + write (iout,*) "Fragment list" + do j=1,nfrag(1) + write (iout,*)"Fragment",j," list",(list_frag(k,j), + & k=1,nlist_frag(j)) + enddo + return + end