PROGRAM DODCA c arbitrary distribution of dipoles for fibril modeling IMPLICIT none integer*4 ns1,ns2,nd,j,i,ii,ia,IG,ix real*8 d1,d2,t1,t2,r,dar,pi,x,y,z,k(3),t,ar,vt(3), 1p,bxo,byo,bzo,bx,by,bz real*8 ,allocatable::v0(:,:),dip(:),v121(:),v122(:),alpha(:), 1ecm(:),vyo(:,:,:),po(:,:) c write(6,6001) 6001 format(/,' DODCA: arbitrary distribution of dipoles', 1 ' for fibril modeling',/,/, 1 ' Input : DODCA.OPT',/, 1 ' Oputput: DC.SC',/) call readopt0(ns1,ns2,nd,d1,d2,t1,t2,r,p) allocate(dip(nd),v121(nd),v122(nd),alpha(nd),ecm(nd)) call readopt1(nd,dip,alpha,ecm,v121,v122) write(6,6002)ns1,ns2,nd,d1,d2,t1,t2,p,r 6002 format(' Current options:',/, 1 ' NS1 number of sites direction 1 ',i12,/, 1 ' NS2 number of sites direction 2 ',i12,/, 1 ' NDIPOLE number of dipoles per site ',i12,/, 1 ' DIST1 d12 distance dir 1, A ',f12.4,/, 1 ' DIST2 d12 distance dir 2, A ',f12.4,/, 1 ' T121 t12 torsion dir 1 (ribbon), deg ',f12.4,/, 1 ' T122 t12 torsion dir 2 (strand), deg ',f12.4,/, 1 ' PITCH lagre helix pitch, % ',f12.4,/, 1 ' RADIUS fold radius, A ',g12.4,/) write(6,6003) 6003 format(' DIPOLE ALPHA ENERGY V121 V122') do 208 i=1,nd 208 write(6,6004)i,dip(i),alpha(i),ecm(i),v121(i),v122(i) 6004 format(i2,g12.4,4f10.2) write(6,*) pi=4.0d0*atan(1.0d0) do 3 i=1,nd 3 alpha(i)=alpha(i)*pi/180.0d0 t1=t1*pi/180.0d0 t2=t2*pi/180.0d0 dar=2.0d0*asin(d1/2.0d0/r) allocate(v0(nd,3),vyo(ns2,nd,3), 1po(ns2,3)) do 1 i=1,nd v0(i,1)=0.0d0 v0(i,2)=dip(i)*sin(alpha(i)) 1 v0(i,3)=dip(i)*cos(alpha(i)) c arbitrary geometry for imaging only: OPEN(2,FILE='TEMP.X') write(2,*)' model geometry of dipoles' write(2,*)ns1*ns2*2 open(3,file='DC.SC') write(3,*)ns1*ns2,' total groups' IG=0 ia=0 c first row: c position of the y-bar: bxo=0.0d0 byo=0.0d0 bzo=0.0d0 do 521 j=1,ns2 IG=IG+1 write(3,*)IG,' . group:' t=t2*(dble(j-1)-dble(ns2-1)/2.0d0) x=0.0d0 y=d2*(dble(j-1)-dble(ns2-1)/2.0d0) z=0.0d0 write(3,1000)x/10.0d0,y/10.0d0,z/10.0d0 1000 format(3f15.6,' nm') c positions in the y-bar: po(j,1)=0.0d0 po(j,2)=y po(j,3)=0.0d0 write(3,*)nd+1,' states' write(3,*) c ground state: write(3,*)'0 0.0' do 201 ii=1,nd write(3,1001)ii,ecm(ii)/8065.54093d0,ecm(ii),1.0d7/ecm(ii) 1001 format(i3,f11.6,' eV =',f11.2,' cm-1 =',f11.2,' nm') c vectors in the bar: vyo(j,ii,1)= cos(t)*v0(ii,1)+sin(t)*v0(ii,3) vyo(j,ii,2)=v0(ii,2) 201 vyo(j,ii,3)=-sin(t)*v0(ii,1)+cos(t)*v0(ii,3) vt(1)=vyo(j,1,1) vt(2)=vyo(j,1,2) vt(3)=vyo(j,1,3) call v2(x,y,z,vt,ia) write(3,*) do 521 ii=1,nd 521 write(3,1002)ii,(vyo(j,ii,ix),ix=1,3),(0.0d0,ix=1,3) 1002 format(i3,6g12.4) c rest: do 2 i=2,ns1 ar=dble(i-1)*dar c bar center of this angle: bx=r*(cos(ar)-1.0d0) by=r*ar*p/100.0d0 bz=r*sin(ar) c shift direction: k(1)=bx-bxo k(2)=by-byo k(3)=bz-bzo call vnorm(k) do 52 j=1,ns2 c old group position vector vs bar center: vt(1)=po(j,1)-bxo vt(2)=po(j,2)-byo vt(3)=po(j,3)-bzo c rotate by t1 around k call rvv(vt,k,t1) x=bx+vt(1) y=by+vt(2) z=bz+vt(3) po(j,1)=x po(j,2)=y po(j,3)=z IG=IG+1 write(3,*)IG,' . group:' write(3,1000)x/10.0d0,y/10.0d0,z/10.0d0 write(3,*)nd+1,' states' write(3,*) c ground state: write(3,*)'0 0.0' do 48 ii=1,nd 48 write(3,1001)ii,ecm(ii)/8065.54093d0,ecm(ii),1.0d7/ecm(ii) write(3,*) do 52 ii=1,nd c previous vector rotated by delta alpha r: vt(1)= cos(dar)*vyo(j,ii,1)+sin(dar)*vyo(j,ii,3) vt(2)=vyo(j,ii,2) vt(3)=-sin(dar)*vyo(j,ii,1)+cos(dar)*vyo(j,ii,3) c vt twisted around k by t1: call rvv(vt,k,t1) c save vyo(j,ii,1)=vt(1) vyo(j,ii,2)=vt(2) vyo(j,ii,3)=vt(3) c write rotated dipole moments: write(3,1002)ii,(vt(ix),ix=1,3),(0.0d0,ix=1,3) c record geom: 52 if(ii.eq.1)call v2(x,y,z,vt,ia) bxo=bx byo=by 2 bzo=bz close(2) close(3) end subroutine vnorm(k) implicit none real*8 k(*),k2,sp k2=dsqrt(sp(k,k)) k(1)=k(1)/k2 k(2)=k(2)/k2 k(3)=k(3)/k2 return end function sp(a,b) implicit none real*8 a(*),b(*),sp sp=a(1)*b(1)+a(2)*b(2)+a(3)*b(3) return end subroutine vp(a,b,c) implicit none real*8 a(*),b(*),c(*) c(1)=a(2)*b(3)-a(3)*b(2) c(2)=a(3)*b(1)-a(1)*b(3) c(3)=a(1)*b(2)-a(2)*b(1) return end subroutine readopt0(ns1,ns2,nd,d1,d2,t1,t2,r,p) implicit none integer*4 ns1,ns2,nd real*8 d1,d2,t1,t2,r,p character*80 s8 ns1=2 ns2=1 nd=1 d1=4.84d0 d2=3.5d0 t1=0.0d0 t2=0.0d0 r=1.0d40 p=1.0d0 open(8,file='DODCA.OPT',status='old') 1 read(8,800,end=99,err=99)s8 800 format(a8) if(s8(1:3).eq.'NS1')read(8,*)ns1 if(s8(1:3).eq.'NS2')read(8,*)ns2 if(s8(1:3).eq.'NDI')read(8,*)nd if(s8(1:5).eq.'DIST1')read(8,*)d1 if(s8(1:5).eq.'DIST2')read(8,*)d2 if(s8(1:4).eq.'T121')read(8,*)t1 if(s8(1:4).eq.'T122')read(8,*)t2 if(s8(1:4).eq.'RADI')read(8,*)r if(s8(1:4).eq.'PITC')read(8,*)p goto 1 99 close(8) return end subroutine readopt1(nd,dip,alpha,e,v121,v122) implicit none integer*4 i,nd real*8 dip(*),alpha(*),v121(*),v122(*),e(*) character*80 s8 do 2 i=1,nd e(i)=1650.0d0 dip(i)=1.0d0 alpha(i)=0.0d0 v121(i)=0.0d0 2 v122(i)=0.0d0 open(8,file='DODCA.OPT',status='old') 1 read(8,800,end=99,err=99)s8 800 format(a8) if(s8(1:3).eq.'DIP')read(8,*)(dip( i),i=1,nd) if(s8(1:3).eq.'ALP')read(8,*)(alpha(i),i=1,nd) if(s8(1:3).eq.'ENE')read(8,*)(e( i),i=1,nd) if(s8(1:4).eq.'V121')read(8,*)(v121(i),i=1,nd) if(s8(1:4).eq.'V122')read(8,*)(v122(i),i=1,nd) goto 1 99 close(8) return end subroutine rvv(v,k,t) c vector v twisted around direction k by t: implicit none real*8 v(*),k(*),sp,t,s,c,kv(3),skv,vx,vy,vz s=sin(t) c=cos(t) call vp(k,v,kv) skv=sp(k,v) vx=v(1) vy=v(2) vz=v(3) v(1)=vx*c+kv(1)*s+k(1)*skv*(1.0d0-c) v(2)=vy*c+kv(2)*s+k(2)*skv*(1.0d0-c) v(3)=vz*c+kv(3)*s+k(3)*skv*(1.0d0-c) return end subroutine v2(x,y,z,v,ia) implicit none real*8 x,y,z,v(*),cx,cy,cz,ox,oy,oz integer*4 ia call vnorm(v) cx=x+v(1)*0.63d0 cy=y+v(2)*0.63d0 cz=z+v(3)*0.63d0 ox=x-v(1)*0.63d0 oy=y-v(2)*0.63d0 oz=z-v(3)*0.63d0 ia=ia+1 write(2,200)6,cx,cy,cz,ia+1 write(2,200)8,ox,oy,oz,ia 200 format(i3,3f18.6,i6,' 0 0 0 0 0 0 0.0') ia=ia+1 return end