program wiggle ! uses perfect energy function to test wiggling use defs; use faces; use Ran_Lux_Mod; implicit none character(len=10) time,date,zone character(len=20) :: rcsbnameofpdb,mypdbinfile,logfile, & mypdboutfile,mydenoutfile,denlogfile character(len=12) :: start character(len=4),dimension(:,:),allocatable :: kindofatom integer(i4b) :: nstart,nres,numberofatoms,nsweeps, & ndenaturing,LWORK,ndensweep, & sweep,length,values(8),intseed,imar,residue, & intervalnum,ndone,zero,LUX,testseed integer(i4b),dimension(:),allocatable :: resnum,CAatomnumber real(dp) :: rmsdistance,den_angle_scale,dtheta, & angle_scale,bestdistance real(dp),dimension(:,:),allocatable :: coord,pdbcoord,coord0 real(dp), dimension(:,:), allocatable :: CA,CA0,pdbCA real(sp), dimension(:), allocatable :: RVec read(5,*); read(5,*) start,nstart,ndone,nres,numberofatoms,nsweeps read(5,*); read(5,*) rcsbnameofpdb,mydenoutfile,mypdbinfile read(5,*); read(5,*) logfile,mypdboutfile,denlogfile read(5,*); read(5,*) ndenaturing,intervalnum,den_angle_scale read(5,*); read(5,*) angle_scale,testseed call date_and_time(date,time,zone,values) if ( nstart == 0 ) then ! normal random start intseed = dot_product(values,scalefactors)*(1+values(8)) else ! stipulated seed intseed = testseed end if zero = 0; LUX = 3 call RLuxGo(LUX,intseed,zero,zero) ! intseed is the seed open(7,file=logfile) write(7,*)'# start, nstart, ndone, nres, numberofatoms, nsweeps' write(7,*)'# ', start, nstart, ndone, nres, numberofatoms, nsweeps write(7,*)'# rcsbnameofpdb, mydenoutfile, mypdbinfile' write(7,*)'# ', rcsbnameofpdb, mydenoutfile, mypdbinfile write(7,*)'# logfile, mypdboutfile, denlogfile' write(7,*)'# ', logfile, mypdboutfile, denlogfile write(7,*)'# ndenaturing, intervalnum, den_angle_scale' write(7,*)'# ', ndenaturing, intervalnum, den_angle_scale write(7,*)'# angle_scale, testseed' write(7,*)'# ', angle_scale, testseed write(7,*)'# WRIGGLING from seed ',intseed,' at ',time,' on ',date close(7) LWORK = 8 allocate(kindofatom(2,numberofatoms),resnum(numberofatoms), & coord(3,numberofatoms),pdbcoord(3,numberofatoms), & CAatomnumber(nres),pdbCA(3,nres),CA(3,nres), & CA0(3,nres),coord0(3,numberofatoms)) call readpdb(rcsbnameofpdb,numberofatoms,kindofatom, & resnum,pdbcoord,CAatomnumber) call CAtrace(nres,pdbcoord,CAatomnumber,pdbCA) if ( start == 'cold' ) then coord = pdbcoord call denature(nres,numberofatoms,kindofatom, & resnum,coord,CAatomnumber, & pdbCA,ndenaturing,ndone, & intervalnum,den_angle_scale,rmsdistance, & ndensweep,denlogfile,mydenoutfile) open(7,file=logfile,position='append') write(7,'(a,i10)') '# LWORK =',LWORK write(7,'(a,f10.6)') '# rms denatured distance =', rmsdistance write(7,'(a,i10)') '# number of denaturing sweeps was',ndensweep close(7) call writepdb(mydenoutfile,numberofatoms,kindofatom, & resnum,coord) stop else if ( start == 'warm' ) then call readpdb(mypdbinfile,numberofatoms,kindofatom, & resnum,coord,CAatomnumber) call denature(nres,numberofatoms,kindofatom, & resnum,coord,CAatomnumber, & pdbCA,ndenaturing,ndone, & intervalnum,den_angle_scale,rmsdistance, & ndensweep,denlogfile,mydenoutfile) open(7,file=logfile,position='append') write(7,'(a,i10)') '# LWORK =',LWORK write(7,'(a,f10.6)') '# rms denatured distance =', rmsdistance write(7,'(a,i10)') '# number of denaturing sweeps was',ndensweep close(7) call writepdb(mydenoutfile,numberofatoms,kindofatom, & resnum,coord) stop else if ( start == 'hot' ) then call readpdb(mypdbinfile,numberofatoms,kindofatom, & resnum,coord,CAatomnumber) end if length = 2*(nres - 2) + 1; allocate(RVec(length)) ! ????? coord0 = coord; call CAtrace(nres,coord,CAatomnumber,CA) CA0 = CA call distance(nres,CA,pdbCA,rmsdistance) bestdistance = rmsdistance do sweep = 1, nsweeps ! start Monte Carlo call RanLux(RVec); imar = 0 ! get vector RVec of random numbers ! we do not rotate about the bonds of CA1, ! which would affect only N1 do residue = 2, nres - 3 imar = imar + 1; dtheta = angle_scale*(RVec(imar) - 0.5_dp) if ( kindofatom(2,CAatomnumber(residue)) == ' PRO' ) then go to 1 ! skip phi; do psi else if ( ( kindofatom(2,CAatomnumber(residue+1)) == ' PRO' ) & .AND. ( kindofatom(2,CAatomnumber(residue+2)) /= ' PRO' ) ) then call phipro(residue,numberofatoms,dtheta,CAatomnumber,coord) else if ( ( kindofatom(2,CAatomnumber(residue+1)) == ' PRO' ) & .AND. ( kindofatom(2,CAatomnumber(residue+2)) == ' PRO' ) ) then call phipropro(residue,numberofatoms,dtheta,CAatomnumber,coord) else call phi_tremble(residue,numberofatoms,dtheta,CAatomnumber,coord) end if call CAtrace(nres,coord,CAatomnumber,CA) call distance(nres,CA,pdbCA,rmsdistance) if ( rmsdistance >= bestdistance ) then coord = coord0; CA = CA0 ! go back to previous state else bestdistance = rmsdistance ! keep new state coord0 = coord; CA0 = CA end if 1 continue imar = imar + 1; dtheta = angle_scale*(RVec(imar) - 0.5_dp ) if ( ( kindofatom(2,CAatomnumber(residue+1)) == ' PRO' ) & .AND. ( kindofatom(2,CAatomnumber(residue+2)) /= ' PRO' ) ) then call pro1psi(residue,numberofatoms,dtheta,CAatomnumber,coord) else if ( ( kindofatom(2,CAatomnumber(residue+1)) /= ' PRO' ) & .AND.( kindofatom(2,CAatomnumber(residue+2)) == ' PRO' ) ) then call pro2psi(residue,numberofatoms,dtheta,CAatomnumber,coord) else if ( ( kindofatom(2,CAatomnumber(residue+1)) == ' PRO' ) & .AND.( kindofatom(2,CAatomnumber(residue+2)) == ' PRO' ) & .AND.( kindofatom(2,CAatomnumber(residue+3)) /= ' PRO' ) ) then call pro12psi(residue,numberofatoms,dtheta,CAatomnumber,coord) else if ( ( kindofatom(2,CAatomnumber(residue+1)) == ' PRO' ) & .AND.( kindofatom(2,CAatomnumber(residue+2)) == ' PRO' ) & .AND.( kindofatom(2,CAatomnumber(residue+3)) == ' PRO' ) ) then call pro123psi(residue,numberofatoms,dtheta,CAatomnumber,coord) else call psi_quiver(residue,numberofatoms,dtheta,CAatomnumber,coord) end if call CAtrace(nres,coord,CAatomnumber,CA) call distance(nres,CA,pdbCA,rmsdistance) if ( rmsdistance >= bestdistance ) then coord = coord0; CA = CA0 ! go back to previous state else bestdistance = rmsdistance ! keep new state coord0 = coord; CA0 = CA end if end do residue = nres - 2 imar = imar + 1; dtheta = angle_scale*(RVec(imar) - 0.5_dp) if ( kindofatom(2,CAatomnumber(residue)) == ' PRO' ) then go to 2 ! skip phi, do psi else if ( ( kindofatom(2,CAatomnumber(residue+1)) == ' PRO' ) & .AND. ( kindofatom(2,CAatomnumber(residue+2)) /= ' PRO' ) ) then call phipro(residue,numberofatoms,dtheta,CAatomnumber,coord) else if ( ( kindofatom(2,CAatomnumber(residue+1)) == ' PRO' ) & .AND. ( kindofatom(2,CAatomnumber(residue+2)) == ' PRO' ) ) then call phipropro(residue,numberofatoms,dtheta,CAatomnumber,coord) else call phi_tremble(residue,numberofatoms,dtheta,CAatomnumber,coord) end if call CAtrace(nres,coord,CAatomnumber,CA) call distance(nres,CA,pdbCA,rmsdistance) if ( rmsdistance >= bestdistance ) then coord = coord0; CA = CA0 ! go back to previous state else bestdistance = rmsdistance ! keep new state coord0 = coord; CA0 = CA end if 2 continue imar = imar + 1; dtheta = angle_scale*(RVec(imar) - 0.5_dp) if ( ( kindofatom(2,CAatomnumber(residue+1)) == ' PRO' ) & .AND. ( kindofatom(2,CAatomnumber(residue+2)) /= ' PRO' ) ) then call pro1psi(residue,numberofatoms,dtheta,CAatomnumber,coord) else if ( ( kindofatom(2,CAatomnumber(residue+1)) /= ' PRO' ) & .AND.( kindofatom(2,CAatomnumber(residue+2)) == ' PRO' ) ) then call pro2psi(residue,numberofatoms,dtheta,CAatomnumber,coord) else if ( ( kindofatom(2,CAatomnumber(residue+1)) == ' PRO' ) & .AND.( kindofatom(2,CAatomnumber(residue+2)) == ' PRO' ) ) then call psi_rot(residue,coord,dtheta,CAatomnumber,numberofatoms) else call psi_quiver(residue,numberofatoms,dtheta,CAatomnumber,coord) end if call CAtrace(nres,coord,CAatomnumber,CA) call distance(nres,CA,pdbCA,rmsdistance) if ( rmsdistance >= bestdistance ) then coord = coord0; CA = CA0 ! go back to previous state else bestdistance = rmsdistance ! keep new state coord0 = coord; CA0 = CA end if residue = nres - 1 ! fix from here on down imar = imar + 1; dtheta = angle_scale*(RVec(imar) - 0.5_dp) if ( kindofatom(2,CAatomnumber(residue)) == ' PRO' ) then go to 3 ! skip phi, do psi else if ( kindofatom(2,CAatomnumber(residue+1)) /= ' PRO' ) then call phi_tremble(residue,numberofatoms,dtheta,CAatomnumber,coord) else call phi_rot(residue,coord,dtheta,Caatomnumber,numberofatoms) end if call CAtrace(nres,coord,CAatomnumber,CA) call distance(nres,CA,pdbCA,rmsdistance) if ( rmsdistance >= bestdistance ) then coord = coord0; CA = CA0 ! go back to previous state else bestdistance = rmsdistance ! keep new state coord0 = coord; CA0 = CA end if 3 continue imar = imar + 1; dtheta = angle_scale*(RVec(imar) - 0.5_dp) call psi_rot(residue,coord,dtheta,CAatomnumber,numberofatoms) call CAtrace(nres,coord,CAatomnumber,CA) call distance(nres,CA,pdbCA,rmsdistance) if ( rmsdistance >= bestdistance ) then coord = coord0; CA = CA0 ! go back to previous state else bestdistance = rmsdistance ! keep new state coord0 = coord; CA0 = CA end if ! we do not rotate about the bonds of CA_nres ! which would only affect the terminal C=O if ( mod(sweep,50) == 0 ) then open(7,file=logfile,position='append') write(7,'(i7,1x,f10.6)') ndone+sweep, bestdistance close(7) call writepdb(mypdboutfile,numberofatoms,kindofatom, & resnum,coord) else if ( nsweeps < 1000 ) then open(7,file=logfile,position='append') write(7,'(i7,1x,f10.6)') ndone+sweep, bestdistance close(7) end if end do ! end of Monte Carlo call writepdb(mypdboutfile,numberofatoms,kindofatom, & resnum,coord) end program wiggle