PROGRAM DODCA IMPLICIT none integer*4 ns1,ns2,ns3,nd,j,i,ii,ia,IG,ix,l,np1,np2,j1,j2 real*8 d1,d2,d3,p3,t1,t2,r,dar,pi,k(3),vt(3),sp,x,y,z,u(3), 1p,bc(3),o(3),di(3),dmi(3),d2i(3),vtt(3),od(3),vy(3),w(3),wo(3), 1l1,l2,bk(3),xl,yl,zl,dj real*8 ,allocatable::v0(:,:),dip(:),v121(:),v122(:),alpha(:,:), 1ecm(:),vyo(:,:),po(:,:),dxyz(:,:),mxyz(:,:),m0(:,:) logical lxyz data vy/0.0d0,1.0d0,0.0d0/ 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,ns3,nd,d1,d2,d3,p3,t1,t2,r,p,np1,np2, 1l1,l2) allocate(dip(nd),v121(nd),v122(nd),alpha(abs(ns2),nd),ecm(nd), 1dxyz(nd,3),mxyz(nd,3)) call readopt1(nd,ns2,dip,alpha,ecm,v121,v122,lxyz,dxyz,mxyz) write(6,6002)ns1,ns2,ns3,np1,np2,nd,l1,l2,d1,d2,d3,t1,t2,p,p3,r 6002 format(' Current options:',/, 1 ' NS1 number of sites direction 1 ',i12,/, 1 ' NS2 number of sites direction 2 ',i12,/, 1 ' NS3 number of sites direction 3 ',i12,/, 1 ' NP1 long range ',i12,/, 1 ' NP2 long range ',i12,/, 1 ' NDIPOLE number of dipoles per site ',i12,/, 1 ' L1 long range A ',F12.4,/, 1 ' L2 long range A ',F12.4,/, 1 ' DIST1 d12 distance dir 1, A ',f12.4,/, 1 ' DIST2 d12 distance dir 2, A ',f12.4,/, 1 ' DIST3 d12 distance dir 3, 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 ' P3 offset in the direction 3, A ',f12.4,/, 1 ' RADIUS fold radius, A ',g12.4,/) if(lxyz)then write(6,6005) 6005 format(' UM DIPOLES x y z ENERGY V121 V122') do 209 i=1,nd 209 write(6,6006)i,(dxyz(i,j),j=1,3),(mxyz(i,j),j=1,3), 1 ecm(i),v121(i),v122(i) 6006 format(i2,6g11.3,3f10.2) else write(6,6003) 6003 format(' DIPOLE ALPHA ENERGY V121 V122') do 208 i=1,nd 208 write(6,6004)i,dip(i),alpha(1,i),ecm(i),v121(i),v122(i) 6004 format(i2,g12.4,4f10.2) endif write(6,*) pi=4.0d0*atan(1.0d0) do 3 i=1,nd do 3 j=1,abs(ns2) 3 alpha(j,i)=alpha(j,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(nd,6),po(abs(ns2),3),m0(nd,3)) c arbitrary geometry for imaging only: OPEN(2,FILE='TEMP.X') write(2,*)' model geometry of dipoles' write(2,*)ns1*abs(ns2)*ns3*2*np1*np2 open(3,file='DC.SC') write(3,*)ns1*abs(ns2)*ns3,' total groups' IG=0 ia=0 c c | z c | c | c | c || di c || o , w c ||------- c --------------------------- y c //\ c // \di-1,k c //od \ / wo c / / c / x c rotate about x for the first point (1,0,0): o=vy c shift i,i-1,i-2: di(1)=0.0d0 di(2)=0.0d0 di(3)=d1 call rvv(di,o,-dar,dmi) call rvv(dmi,o,-dar,d2i) c bar direction w(1)=0.0d0 w(2)=1.0d0 w(3)=0.0d0 call rvv(w,di,-t1,wo) c position of bar center: bc(1)=0.0d0 bc(2)=0.0d0 bc(3)=-d1 c do 2 i=1,ns1 c c momentary coordinate system: c o - rotation vector: call vp(d2i,dmi,o) if(sp(o,o).lt.1.0d-30)o=vy call vnorm(o) c shift direction: k=dmi call vnorm(k) c perpendicular direction: call vp(di,o,od) call vnorm(od) c di: new shift from di-1: call rvv(dmi,o,dar,di) c conserve length: call rnorm(di,d1) c apply pitch: do 53 ix=1,3 53 di(ix)=di(ix)+p/100.0d0*d1*o(ix) c bc: bar center of this angle: do 102 ix=1,3 102 bc(ix)=bc(ix)+di(ix) c update bar direction: k rotation by t1 and o rotation by dar call rvv(wo,k,t1,vt) call rvv(vt,o,dar,w) call vnorm(w) c vector completing the bar local system w u bk: call vp(w,di,u) call vnorm(u) call vp(u,w,bk) do 52 j=1,abs(ns2) c zero direction vectors in w u bk coordinates (bkw plane): do 1 ii=1,nd do 1 ix=1,3 if(lxyz)then v0(ii,ix)=dxyz(ii,2)*w(ix)+dxyz(ii,3)*k(ix)+dxyz(ii,1)*u(ix) m0(ii,ix)=mxyz(ii,2)*w(ix)+mxyz(ii,3)*k(ix)+mxyz(ii,1)*u(ix) else v0(ii,ix)=dip(ii)*(sin(alpha(j,ii))*w(ix)+cos(alpha(j,ii))*k(ix)) m0(ii,ix)=0.0d0 endif 1 continue c group positions in the bar along w: do 104 ix=1,3 104 po(j,ix)=bc(ix)+d2*(dble(j-1)-dble(abs(ns2)-1)/2.0d0)*w(ix) c new vectors = zero v0 rotated by t around w: if(mod(ns2,2).eq.0)then c even number,zero for j = ns2/2 dj=dble(j-abs(ns2)/2) else c odd number,zero for j = (ns2+1)/2 dj=dble(j-1-(abs(ns2)-1)/2) endif do 106 ii=1,nd do 1061 ix=1,3 1061 vt(ix)=v0(ii,ix) call rvv(vt,w,t2*dj,vtt) do 1063 ix=1,3 1063 vyo(ii,ix)=vtt(ix) do 1062 ix=1,3 1062 vt(ix)=m0(ii,ix) call rvv(vt,w,t2*dj,vtt) do 106 ix=1,3 106 vyo(ii,3+ix)=vtt(ix) do 52 l=1,ns3 xl=po(j,1)+dble(l-1)*(p3*w(1)+d3*u(1)) yl=po(j,2)+dble(l-1)*(p3*w(2)+d3*u(2)) zl=po(j,3)+dble(l-1)*(p3*w(3)+d3*u(3)) do 52 j1=1,np1 do 52 j2=1,np2 IG=IG+1 c positions: just shift in the w and u direction: x=xl+dble(j1-1)*l1*w(1)+dble(j2-1)*l2*u(1) y=yl+dble(j1-1)*l1*w(2)+dble(j2-1)*l2*u(2) z=zl+dble(j1-1)*l1*w(3)+dble(j2-1)*l2*u(3) write(3,*)IG,' . group:' write(3,1000)x/10.0d0,y/10.0d0,z/10.0d0 1000 format(3f15.6,' nm') 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) 1001 format(i3,f11.6,' eV =',f11.2,' cm-1 =',f11.2,' nm') write(3,*) c vectors: the same for all l j1 j2: do 52 ii=1,nd do 107 ix=1,3 107 vt(ix)=vyo(ii,ix) c record geom for ii=1 and this dipole: if(ii.eq.1)call v2(x,y,z,vt,ia) c write rotated dipole moments: 52 write(3,1002)ii,(vt(ix),ix=1,3),(vyo(ii,3+ix),ix=1,3) 1002 format(i3,6g12.4) wo=w d2i=dmi 2 dmi=di 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 subroutine rnorm(k,o) implicit none real*8 k(*),k2,sp,o k2=dsqrt(sp(k,k)) k(1)=k(1)/k2*o k(2)=k(2)/k2*o k(3)=k(3)/k2*o 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,ns3,nd,d1,d2,d3,p3,t1,t2,r,p, 1np1,np2,l1,l2) implicit none integer*4 ns1,ns2,ns3,nd,np1,np2 real*8 d1,d2,d3,p3,t1,t2,r,p,l1,l2 character*80 s8 ns1=2 ns2=1 ns3=1 np1=1 np2=1 nd=1 d1=4.84d0 d2=3.5d0 d3=11.0d0 p3=1.5d0 t1=0.0d0 t2=0.0d0 r=1.0d40 p=1.0d0 l1=100 l2=100 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.'NS3')read(8,*)ns3 if(s8(1:3).eq.'NP1')read(8,*)np1 if(s8(1:3).eq.'NP2')read(8,*)np2 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:5).eq.'DIST3')read(8,*)d3 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 if(s8(1:2).eq.'P3')read(8,*)p3 if(s8(1:2).eq.'L1')read(8,*)l1 if(s8(1:2).eq.'L2')read(8,*)l2 goto 1 99 close(8) return end subroutine readopt1(nd,ns2,dip,alpha,e,v121,v122, 1lxyz,dxyz,mxyz) implicit none integer*4 i,j,ns2,nd real*8 dip(*),alpha(ns2,nd),v121(*),v122(*),e(*), 1dxyz(nd,3),mxyz(nd,3) logical lxyz character*80 s8 do 2 i=1,nd e(i)=1650.0d0 dip(i)=1.0d0 do 22 j=1,ns2 22 alpha(j,i)=0.0d0 v121(i)=0.0d0 2 v122(i)=0.0d0 lxyz=.false. 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) c dipoles in cartesian coordinates: if(s8(1:3).eq.'DXY')then do 5 i=1,nd 5 read(8,*)(dxyz(i,j),j=1,3),(mxyz(i,j),j=1,3) lxyz=.true. endif if(s8(1:3).eq.'ALP')then if(ns2.lt.0)then do 4 j=1,abs(ns2) 4 read(8,*)(alpha(j,i),i=1,nd) else read(8,*)(alpha(1,i),i=1,nd) do 3 j=2,ns2 do 3 i=1,nd 3 alpha(j,i)=alpha(1,i) endif endif 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,vn) c vector v twisted around direction k by t: implicit none real*8 v(*),k(*),sp,t,s,c,kv(3),skv,vn(*) s=sin(t) c=cos(t) call vp(k,v,kv) skv=sp(k,v) vn(1)=v(1)*c+kv(1)*s+k(1)*skv*(1.0d0-c) vn(2)=v(2)*c+kv(2)*s+k(2)*skv*(1.0d0-c) vn(3)=v(3)*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(3),ox,oy,oz,u(3),d integer*4 ia u=v d=1.26d0 call vnorm(u) ox=x-u(1)*d oy=y-u(2)*d oz=z-u(3)*d ia=ia+1 write(2,200)6, x, y, z,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