반응형
MATLAB
Code
xyz = readdcd('C.dcd' , 1: 303);
dcd 파일에서 프레임 번호마다 첫 번째 원자부터 마지막 303번째 원자까지의 x좌표, y좌표, z좌표를 xyz에 저장한다.
atomlist = 5:10:305;
xyz2 = readdcd('C.dcd,',atomlist);
dcd 파일에서 프레임 번호마다 첫번째 alpha-Carbon 부터 마지막 alpha Carbon까지의 x좌표, y좌표, z좌표를 xyz2에 저장한다.
alanin_x = xyz2(:, 1:3:end);
alanin_y = xyz2(:, 2:3:end);
alanin_z = xyz2(:, 3:3:end);
xyz2에서 alpha Carbon 들의 x좌표, y좌표, z좌표를 각각 alanin_x, alanin_y, alanin_x에 저장한다.
%a=1;
%while a<=30
% while b<=3
% centerX=centerX+alanin_x(a,b);
% centerY=centerY+alanin_x(a,b);
% centerZ=centerZ+alanin_x(a,b);
% end
% centerX=centerX/3;
% centerY=centerX/3;
% centerZ=centerX/3;
% a=a+1;
%end
각 슬라이스의 무게중심의 x 좌표, y좌표, z좌표를 각각 centerX, centerY, centerZ에 저장할려고 했으나 프레임 번호를 고려하지 않아 오류가 생겨 밑에 새로 썼다.
frame=5000;
slice=7;
centerX=zeros(frame,slice);
centerY=zeros(frame,slice);
centerZ=zeros(frame,slice);
a=1;
b=1;
center=zeros(frame,slice,3);
while a<=frame
while b<=slice
%if rem(b,4)~=0
centerX(a,b)=alanin_x(a,b*4-3)+alanin_x(a,b*4-2)+alanin_x(a,b*4-1);
centerY(a,b)=alanin_y(a,b*4-3)+alanin_y(a,b*4-2)+alanin_y(a,b*4-1);
centerZ(a,b)=alanin_z(a,b*4-3)+alanin_z(a,b*4-2)+alanin_z(a,b*4-1);
%end
centerX(a,b)=centerX(a,b)/3;
centerY(a,b)=centerY(a,b)/3;
centerZ(a,b)=centerZ(a,b)/3;
프레임 5000까지 각 프레임에 대한 각 슬라이스의 각 꼭지점의 x좌표,y좌표,z좌표 값의 합을 구한다.
center(a,b,:)=[centerX(a,b),centerY(a,b),centerZ(a,b)];
각 슬라이스의 X, Y, Z좌표 각각의 합을 각각 3으로 나누어 X, Y, Z좌표의 질량 중심을 구한 뒤 하나의 이차원 행렬에 저장한다.
b=b+1;
end
a=a+1;
b=1;
end
a=1;
b=1;
centerVecterX=zeros(frame,slice-1);
centerVecterY=zeros(frame,slice-1);
centerVecterZ=zeros(frame,slice-1);
centerDistance=zeros(frame,slice-1);
centerVectorTotal=zeros(frame,slice-1,3);
while a<=frame
while b<=slice-1
centerVectorX(a,b)=centerX(a,b+1)-centerX(a,b);
centerVectorY(a,b)=centerY(a,b+1)-centerY(a,b);
centerVectorZ(a,b)=centerZ(a,b+1)-centerZ(a,b);
centerVectorTotal(a,b,:)=[centerVectorX(a,b),centerVectorY(a,b),centerVectorZ(a,b)];
centerDistance(a,b)=norm(centerVectorTotal(a,b));
b=b+1;
end
각각의 X, Y, Z 좌표의 센터에서 그 이전의 X, Y, Z좌표의 센터를 빼어 각각의 센터 벡터를 구한다. 그 다음 각각의 센터벡터를 하나의 이차원 행렬에 저장한 뒤 norm 시켜 각각의 X, Y, Z좌표의 센터벡터를 단위가 1인 벡터로 변환한다.
b=1;
a=a+1;
end
deltaS=zeros(frame,1);
distanceSum=0;
for m=1:frame
for n=1:(slice-2)
distanceSum=distanceSum+centerDistance(m,n);
end
deltaS(m,1)=distanceSum/(slice-1);
distanceSum=0;
end
%deltaS=mean(deltaS);
sliceVector=zeros(frame,slice-1,3);
e3=zeros(frame,slice-1,3);
for m=1:frame
for n=1:(slice-1)
sliceVector(m,n,:)=center(m,n+1,:)-center(m,n,:);
%e3(m,n,:)=sliceVector(m,n,:)/norm(sliceVector(m,n,:));
e3(m,n,:)=sliceVector(m,n,:)/norm(sliceVector(m,n));
end
end
ePrimeX=zeros(frame,slice);
ePrimeY=zeros(frame,slice);
ePrimeZ=zeros(frame,slice);
ePrimeVectorX=zeros(frame,slice);
ePrimeVectorY=zeros(frame,slice);
ePrimeVectorZ=zeros(frame,slice);
ePrime=zeros(frame,slice,3);
ePrimeVectorTotal=zeros(frame,slice,3);
for m=1:frame
for n=1:slice
a=4*n-3;
ePrimeX(m,n)=(alanin_x(m,a)-alanin_x(m,a+2))/2;
ePrimeY(m,n)=(alanin_y(m,a)-alanin_y(m,a+2))/2;
ePrimeZ(m,n)=(alanin_z(m,a)-alanin_z(m,a+2))/2;
ePrimeVectorX(m,n)=ePrimeX(m,n)-centerX(m,n,1);
ePrimeVectorY(m,n)=ePrimeY(m,n)-centerY(m,n,1);
ePrimeVectorZ(m,n)=ePrimeZ(m,n)-centerZ(m,n,1);
ePrime(m,n,:)=[ePrimeX(m,n),ePrimeY(m,n),ePrimeZ(m,n)];
ePrimeVectorTotal(m,n,:)=[ePrimeVectorX(m,n),ePrimeVectorY(m,n),ePrimeVectorZ(m,n)];
end
end
프레임에 다하여 각 슬라이스의 꼭지점을 이루는 알파카본 중 1번과 3번의 중심을 구한다.
%ePrimeTotalX=zeros(frame,slice);
%ePrimeTotalY=zeros(frame,slice);
%ePrimeTotalZ=zeros(frame,slice);
%a=1;
%b=1;
%while a<=frame
% while b<=slice
% end
%end
e1Vector=zeros(frame,slice-1,3);
e2Vector=zeros(frame,slice-1,3);
e3Vector=zeros(frame,slice-1,3);
e1VectorCross=zeros(frame,slice-1,3);
e2VectorCross=zeros(frame,slice-1,3);
e1Vector0=zeros(frame,slice-1,3);
e2Vector0=zeros(frame,slice-1,3);
e3Vector0=zeros(frame,slice-1,3);
for m=1:frame
for n=1:(slice-1)
%CenterVectorTotal
%ePrimeVectorTotal
%e2VectorCross(m,n,:)=cross(centerVectorTotal(m,n,:), ePrimeVectorTotal(m,n,:));
%e1VectorCross(m,n,:)=cross(e2VectorCross(m,n,:),centerVectorTotal(m,n,:));
%{
e2VectorCross(m,n)=cross(centerVectorTotal(m,n), ePrimeVectorTotal(m,n));
e1VectorCross(m,n)=cross(e2VectorCross(m,n),centerVectorTotal(m,n));
e2Vector0(m,n)=e2VectorCross(m,n)/norm(e2VectorCross(m,n));
e1Vector0(m,n)=e1VectorCross(m,n)/norm(e1VectorCross(m,n));
e3Vector0(m,n)=centerVectorTotal(m,n)/norm(centerVectorTotal(m,n))
e1Vector(m,n,1)=e1Vector0(m,n,1);
e1Vector(m,n,2)=e1Vector0(m,n,2);
e1Vector(m,n,3)=e1Vector0(m,n,3);
e2Vector(m,n,1)=e2Vector0(m,n,1);
e2Vector(m,n,2)=e2Vector0(m,n,2);
e2Vector(m,n,3)=e2Vector0(m,n,3);
e3Vector(m,n,1)=e3Vector0(m,n,1);
e3Vector(m,n,2)=e3Vector0(m,n,2);
e3Vector(m,n,3)=e3Vector0(m,n,3);
%}
앞에서 구한 e1VectorCross, e2VectorCross,e3VectorCross를 단위 벡터 처리 해준다.
centerVectorTotal0=[centerVectorX(m,n),centerVectorY(m,n),centerVectorZ(m,n)];
ePrimeVectorTotal0=[ePrimeVectorX(m,n),ePrimeVectorY(m,n),ePrimeVectorZ(m,n)];
e2VectorCrossAns=cross(centerVectorTotal0, ePrimeVectorTotal0);
e1VectorCrossAns=cross(e2VectorCrossAns, centerVectorTotal0);
%크기가 1인 벡터로 조정
e1VectorNorm=e1VectorCrossAns/norm(e1VectorCrossAns);
e2VectorNorm=e2VectorCrossAns/norm(e2VectorCrossAns);
e3VectorNorm=centerVectorTotal0/norm(centerVectorTotal0);
앞에서 구한 CenterVectorTotal, e1VectorCrossAns, e2VectorCrossAns를 크기가 1인 단위벡터로 처리한다.
e1Vector(m,n,1)=e1VectorNorm(1);
e1Vector(m,n,2)=e1VectorNorm(2);
e1Vector(m,n,3)=e1VectorNorm(3);
e2Vector(m,n,1)=e2VectorNorm(1);
e2Vector(m,n,2)=e2VectorNorm(2);
e2Vector(m,n,3)=e2VectorNorm(3);
e3Vector(m,n,1)=e3VectorNorm(1);
e3Vector(m,n,2)=e3VectorNorm(2);
e3Vector(m,n,3)=e3VectorNorm(3);
end
end
단위 벡터 처리한 벡터에서 x,y,z의 값을 추출하여 따로 저장한다
U=zeros(frame,slice-2,3,3);
for m=1:frame
for n=1:(slice-2)
%c=1
%for c<=3
% d=1
% for d<=3
% U(a,b,c,d)=dot(ecVector(a,b+1,:),edVector(a,b,:));
% end
%end
%변수이름에 변수사용 저장법 O 변수를 불러올때 변수 사용?
U(m,n,1,1)=dot(e1Vector(m,n+1,:),e1Vector(m,n,:));
U(m,n,2,1)=dot(e2Vector(m,n+1,:),e1Vector(m,n,:));
U(m,n,3,1)=dot(e3Vector(m,n+1,:),e1Vector(m,n,:));
U(m,n,1,2)=dot(e1Vector(m,n+1,:),e2Vector(m,n,:));
U(m,n,2,2)=dot(e2Vector(m,n+1,:),e2Vector(m,n,:));
U(m,n,3,2)=dot(e3Vector(m,n+1,:),e2Vector(m,n,:));
U(m,n,1,3)=dot(e1Vector(m,n+1,:),e3Vector(m,n,:));
U(m,n,2,3)=dot(e2Vector(m,n+1,:),e3Vector(m,n,:));
U(m,n,3,3)=dot(e3Vector(m,n+1,:),e3Vector(m,n,:));
end
end
e1, e2, e3 벡터 세가지 중 (한개 벡터 중복 내적도 계산) 벡터를 2개씩 묶어서 선택한 2개의 벡터끼리 각각 내적하여 U값을 구한다.
%UfrmaeOne=zeros(1,,,)
UframeOne=U(1,:,:,:);
smallU=zeros(3,3);
smallOmega=zeros(frame,slice-2,3,3);
for m=1:frame
for n=1:(slice-2)
smallU(:,:)=U(m,n,:,:);
%smallOmega(m,n,:,:)=(-logm(abs(smallU))/deltaS(m,1));
%In logm (line 50)
%In son (line 198)
smallOmega(m,n,:,:)=(-logm(smallU))/deltaS(m,1);
%smallOmega(m,n,:,:)=smallU;
end
end
이때 U=e^오메가*델타S라는 식에 따라 각 프레임에 대한 오메가 값을 구하기 위하여 앞에서 구한 U를 로그 처리한 후 델타S로 나누어준다.
omega1=zeros(frame,slice-2);
omega2=zeros(frame,slice-2);
omega3=zeros(frame,slice-2);
for m=1:frame
for n=1:(slice-2)
omega1(m,n)=smallOmega(m,n,2,3);
omega2(m,n)=smallOmega(m,n,3,1);
omega3(m,n)=smallOmega(m,n,1,2);
end
end
우리에게 필요한 값인 (2,3), (3,1),(1,2)번째 오메가 값을 추출한다.
%subplot(2,2,1)
%histogram(omega1)
%subplot(2,2,2)
%histogram(omega2)
%subplot(2,2,3)
%histogram(omega3)
%plot3(center(:,:,1))
%plot3(center(:,:,2))
%plot3(center(:,:,3))
%
% for x = 1:10
% disp(x)
% end
subplot(4,3,1)
histogram(omega1(:,1))
subplot(4,3,2)
histogram(omega2(:,1))
subplot(4,3,3)
histogram(omega3(:,1))
%정규화작업과 alpha값 구하기
막대의 X간격을 설정한다. 이때 미리 히스토그램으로 본 omega의 수치들을 확인 한 후 정한다.(omega1, omega2의 경우 0의 중심으로) 좌우대칭이 된다. omega3의 경우 helix가 나선구조이므로 0이 아닌 -값을 중심으로 좌우대칭이 된다.
omega1DeltaLeft=-1.00;
omega1DeltaRight=1.00;
omega2DeltaLeft=-1.00;
omega2DeltaRight=1.00;
omega3DeltaLeft=-1.50;
omega3DeltaRight=0.50;
deltaX1Size=(omega1DeltaRight-omega1DeltaLeft)/2000;
deltaX2Size=(omega2DeltaRight-omega2DeltaLeft)/2000;
deltaX3Size=(omega3DeltaRight-omega3DeltaLeft)/2000;
count1=zeros(slice-2,2000);
count2=zeros(slice-2,2000);
count3=zeros(slice-2,2000);
각 슬라이스별로 오메가 값 붙포를 히스토그램으로 나타낸다.(히스토그램으로 나타낸다는 것을 count변수를 활용한다. )
for m=1:frame
for n=1:slice-2
for c=1:2000
if ((omega1DeltaLeft+c*deltaX1Size)<=omega1(m,n)) && (omega1(m,n)<(omega1DeltaLeft+(c+1)*deltaX1Size))
count1(n,c)=count1(n,c)+1;
end
if ((omega2DeltaLeft+c*deltaX2Size)<=omega2(m,n)) && (omega2(m,n)<(omega2DeltaLeft+(c+1)*deltaX2Size))
count2(n,c)=count2(n,c)+1;
end
if ((omega3DeltaLeft+c*deltaX3Size)<=omega3(m,n)) && (omega3(m,n)<(omega3DeltaLeft+(c+1)*deltaX3Size))
count3(n,c)=count3(n,c)+1;
end
end
end
end
maxcount1=zeros(slice-2,1);
maxcount2=zeros(slice-2,1);
maxcount3=zeros(slice-2,1);
takecount1=zeros(slice-2,2000);
takecount2=zeros(slice-2,2000);
takecount3=zeros(slice-2,2000);
for n=1:slice-2
for c=1:2000
maxcount1=max(count1(n,:));
maxcount2=max(count2(n,:));
maxcount3=max(count3(n,:));
if count1(n,c)>=maxcount1(n,c)/10
takecount1(n,c)=takecount1(n,c)+1;
end
if count2(n,c)>=maxcount2(n,c)/10
takecount2(n,c)=takecount2(n,c)+1;
end
if count3(n,c)>=maxcount3(n,c)/10
takecount3(n,c)=takecount3(n,c)+1;
end
end
end
sum1=zeros(slice-2,1);
sum2=zeros(slice-2,1);
sum3=zeros(slice-2,1);
temp1=0;
temp2=0;
temp3=0;
for n=1:slice-2
for o=1:2000ㄱ
temp1=count1(n,o)*deltaX1Size;
temp2=count2(n,o)*deltaX2Size;
temp3=count3(n,o)*deltaX3Size;
sum1(n,1)=sum1(n,1)+temp1;
sum2(n,1)=sum2(n,1)+temp2;
sum3(n,1)=sum3(n,1)+temp3;
end
end
count1Norm=zeros(slice,2000);
count2Norm=zeros(slice,2000);
count3Norm=zeros(slice,2000);
for m=1:slice-2
count1Norm(m,:)=count1(m,:)/sum1(m,1);
count2Norm(m,:)=count2(m,:)/sum2(m,1);
count3Norm(m,:)=count3(m,:)/sum3(m,1);
end
%{
subplot(4,3,4)
plot(count1Norm)
subplot(4,3,5)
plot(count2Norm)
subplot(4,3,6)
plot(count3Norm)
subplot(4,3,7)
histogram(count1Norm)
subplot(4,3,8)
histogram(count2Norm)
subplot(4,3,9)
histogram(count3Norm)
%}
위에서 구한 값들을 figure 화면 위에 그래프를 그린다.반응형
댓글