X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fmap.F90;fp=source%2Funres%2Fmap.F90;h=b91d43e105fbe212d23ab11be61e7deb91ef301e;hb=2611e37f82b576a1366c2d78fce87c1a55852037;hp=0000000000000000000000000000000000000000;hpb=a1e9d5a6b704c90ebd10f86a655780212a880dce;p=unres4.git diff --git a/source/unres/map.F90 b/source/unres/map.F90 new file mode 100644 index 0000000..b91d43e --- /dev/null +++ b/source/unres/map.F90 @@ -0,0 +1,191 @@ + module map_ +!----------------------------------------------------------------------------- + use io_units + use names + use geometry_data + use control_data + use energy_data + use map_data + implicit none +!----------------------------------------------------------------------------- +! +! +!----------------------------------------------------------------------------- + contains +!----------------------------------------------------------------------------- +! check_sc_map.f +!----------------------------------------------------------------------------- + subroutine check_sc_map +! Subroutine is checking if the fitted function which describs sc_rot_pot +! is correct, printing, alpha,beta, energy, data - for some known theta. +! theta angle is read from the input file. Sc_rot_pot are printed +! for the second residue in sequance. + use geometry, only:chainbuild + use energy, only:vec_and_deriv,esc +! include 'DIMENSIONS' +! include 'COMMON.VAR' +! include 'COMMON.GEO' +! include 'COMMON.INTERACT' + real(kind=8) :: xx,yy,zz,al,om + real(kind=8) :: escloc, escloc_min + real(kind=8),dimension(50000) :: escloc_ene,alph_plot,beta_plot + integer,dimension(5000) :: al_plot,be_plot + integer :: iialph, iibet,it + +!el local variables + integer :: i,j + + write (2,*) "Side-chain-rotamer potential energy map!!!!" + escloc_min = 1000000.00 +! it=itype(2) + i = 0 + do iialph=0,18 + do iibet=-18,18 + i = i + 1 + al = iialph*10.0d0*deg2rad + om = iibet*10.0d0*deg2rad + zz = dcos(al) + xx = -dsin(al)*dcos(om) + yy = -dsin(al)*dsin(om) + alph(2)=dacos(xx) + omeg(2)=-datan2(zz,yy) + al_plot(i)=alph(2)*rad2deg + be_plot(i)=omeg(2)*rad2deg +! write(2,*) alph(2)*rad2deg, omeg(2)*rad2deg + alph_plot(i) = al*rad2deg + beta_plot(i) = om*rad2deg + call chainbuild + call vec_and_deriv + call esc(escloc) + escloc_ene(i) = escloc + if (escloc_min.gt.escloc_ene(i)) escloc_min=escloc_ene(i) + enddo + enddo +! write (2,*) "escloc_min = ", escloc_min + print *,"i",i + do j = 1,i + write (2,'(3f10.3,2i9,f12.5)') alph_plot(j), & + beta_plot(j),theta(3)*rad2deg, al_plot(j),be_plot(j),& + escloc_ene(j) !- escloc_min + enddo + return + end subroutine check_sc_map +!----------------------------------------------------------------------------- +! map.f +!----------------------------------------------------------------------------- + subroutine map + + use geometry, only:chainbuild,geom_to_var + use energy + use minimm, only:minimize +! implicit real*8 (a-h,o-z) +! include 'DIMENSIONS' +! include 'COMMON.MAP' +! include 'COMMON.VAR' +! include 'COMMON.GEO' +! include 'COMMON.DERIV' +! include 'COMMON.IOUNITS' +! include 'COMMON.NAMES' +! include 'COMMON.CONTROL' +! include 'COMMON.TORCNSTR' + real(kind=8) :: energia(0:n_ene) + character(len=5) :: angid(4)=reshape((/'PHI ','THETA','ALPHA',& + 'OMEGA'/),shape(angid)) + real(kind=8) :: ang_list(10) +!el real(kind=8),dimension(6*nres) :: g,x !(maxvar) (maxvar=6*maxres) + real(kind=8),dimension(:),allocatable :: g,x !(maxvar) (maxvar=6*maxres) + integer :: nn(10) +!el local variables + integer :: i,iii,ii,j,k,nmax,ntot,nf,iretcode,nfun + real(kind=8) :: etot,gnorm!,fdum + integer,dimension(1) :: uiparm + real(kind=8),dimension(1) :: urparm + allocate(x(6*nres),g(6*nres)) + + write (iout,'(a,i3,a)')'Energy map constructed in the following ',& + nmap,' groups of variables:' + do i=1,nmap + write (iout,'(2a,i3,a,i3)') angid(kang(i)),' of residues ',& + res1(i),' to ',res2(i) + enddo + nmax=nstep(1) + do i=2,nmap + if (nmax.lt.nstep(i)) nmax=nstep(i) + enddo + ntot=nmax**nmap + iii=0 + write (istat,'(1h#,a14,29a15)') (" ",k=1,nmap),& + (ename(print_order(k)),k=1,nprint_ene),"ETOT","GNORM" + do i=0,ntot-1 + ii=i + do j=1,nmap + nn(j)=mod(ii,nmax)+1 + ii=ii/nmax + enddo + do j=1,nmap + if (nn(j).gt.nstep(j)) goto 10 + enddo + iii=iii+1 +!d write (iout,*) i,iii,(nn(j),j=1,nmap) + do j=1,nmap + ang_list(j)=ang_from(j) & + +(nn(j)-1)*(ang_to(j)-ang_from(j))/nstep(j) + do k=res1(j),res2(j) + goto (1,2,3,4), kang(j) + 1 phi(k)=deg2rad*ang_list(j) + if (minim) phi0(k-res1(j)+1)=deg2rad*ang_list(j) + goto 5 + 2 theta(k)=deg2rad*ang_list(j) + goto 5 + 3 alph(k)=deg2rad*ang_list(j) + goto 5 + 4 omeg(k)=deg2rad*ang_list(j) + 5 continue + enddo ! k + enddo ! j + call chainbuild + if (minim) then + call geom_to_var(nvar,x) + call minimize(etot,x,iretcode,nfun) + print *,'SUMSL return code is',iretcode,' eval ',nfun +! call intout + else + call zerograd + call geom_to_var(nvar,x) + endif + call etotal(energia) + etot = energia(0) + nf=1 + nfl=3 + call gradient(nvar,x,nf,g,uiparm,urparm,fdum) + gnorm=0.0d0 + do k=1,nvar + gnorm=gnorm+g(k)**2 + enddo + etot=energia(0) + + gnorm=dsqrt(gnorm) +! write (iout,'(6(1pe15.5))') (ang_list(k),k=1,nmap),etot,gnorm + write (istat,'(30e15.5)') (ang_list(k),k=1,nmap),& + (energia(print_order(ii)),ii=1,nprint_ene),etot,gnorm +! write (iout,*) 'POINT',I,' ANGLES:',(ang_list(k),k=1,nmap) +! call intout +! call enerprint(energia) + 10 continue + enddo ! i +! deallocate(x,g) + return + end subroutine map +!----------------------------------------------------------------------------- + subroutine alloc_map_arrays + +! commom.map +! common /mapp/ + allocate(kang(nmap),res1(nmap),res2(nmap),nstep(nmap)) !(maxvar) + allocate(ang_from(nmap),ang_to(nmap)) !(maxvar) + + return + end subroutine alloc_map_arrays +!----------------------------------------------------------------------------- +!----------------------------------------------------------------------------- + end module map_