FORTRAN 90 and Python KRISS http://ihlee.kriss.re.kr/compphys/f90module.htm http://ihlee.kriss.re.kr/compphys/pfphysics.htm Stay hungry, stay foolish.
Language is the house of the truth of Being. -Martin Heidegger.
Computer programming (often shortened to programming or coding) is the process of writing, testing, debugging/troubleshooting, and maintaining the source code of computer programs. A programming language is an artificial language that can be used to control the behavior of a machine, particularly a computer. A prominent purpose of programming languages is to provide instructions to a computer. Efficiency, reliability, robustness, portability, readability :!,. Fortran 90 and Python. -Martin Heidegger 1.. 2.. 3.. debugging Science..
Debugging is a very important task for every programmer, because an erroneous program is often useless. Bugs - Compile-time Errors Bugs - Run-time Errors 1945 9 9 bug; `debugging`..
Why FORTRAN 90? FORTRAN 90 1.. ( ;.) 2.. (, 113, ) 3.. (,, ) 4.. (,, ) 5.. ( 77.) ()?: (object), ( /...)... 90.... 1. data abstraction -- user-defined types 2. data hiding -- private and public attributes 3. encapsulation -- modules and data hiding facilities 4. inheritance and extensibility -- generic procedures 5. polymorphism -- generic overloading 6. reusability -- modules, ( ) John Backus
!234567890 program HW implicit none write(6,*) 'Hello World!' 77. C c *! Comment F90.! comment. end program hw ihlee@cetus0:~/lecture_programming$ a.out Hello World! ihlee@cetus0:~/lecture_programming$ 90 : http://incredible.egloos.com/3044862
. FORTRNA 90.
abc=1.0d0 abc 1.d0. 1.e0 (single precision) vs 1.d0 (double precision) 1.d0. (,,,,,,.) () Implicit Typing Undeclared variables have an implicit type, if first letter is I, J, K, L, M or N then type is INTEGER; any other letter then type is REAL s. Operator precedence pi*radius**2 = pi*(radius**2) x*y/z+a= (x*y/z)+a
, Intrinsic Logical Operations A LOGICAL or boolean expression returns a.true. /.FALSE. result. The following are valid with LOGICAL operands,.not. --.TRUE. if operand is.false...and. --.TRUE. if both operands are.true.;.or. --.TRUE. if at least one operand is.true.;.eqv. --.TRUE. if both operands are the same;.neqv. --.TRUE. if both operands are different. For example, if T is.true. and F is.false. FORTRAN 77 do while.. while (logical expr) do statements enddo do while (logical expr) statements enddo integer n n = 1 10 if (n.le. 100) then n = 2*n write (*,*) n goto 10 endif
IF (I > 17) THEN Print*, "I > 17" ELSE Print*, "I <= 17" END IF IF (I > 17) THEN Print*, "I > 17" ELSEIF (I == 17) THEN Print*, "I == 17" ELSE Print*, "I < 17" END IF outa: IF (a.ne. 0) THEN PRINT*, "a /= 0" IF (c.ne. 0) THEN PRINT*, "a /= 0 AND c /= 0" ELSE PRINT*, "a /= 0 BUT c == 0" ENDIF ELSEIF (a.gt. 0) THEN outa PRINT*, "a > 0" ELSE outa PRINT*, "a must be < 0" ENDIF outa [name:] IF (expression) THEN block ENDIF [name] [name:] IF (expression) THEN block1 ELSE [name] block2 ENDIF [name] elseif.. IF ( n == 1 ) THEN discount = 0.0d0 ELSEIF ( n <=5 ) THEN discount = 0.1d0 ELSEIF ( n>5.and. n <=10) THEN discount = 0.15d0 ELSE discount = 0.25d0 ENDIF
INTEGER :: month season: SELECT CASE( month ) CASE(4,5) WRITE(*,*) `Spring' CASE(6,7) WRITE(*,*) `Summer' CASE(8:10) WRITE(*,*) `Autumn' CASE(11,1:3,12) WRITE(*,*) `Winter' CASE DEFAULT WRITE(*,*) `not a month' END SELCET season Conditional statements IF statement IF construct CASE construct Repetition DO constructs GOTO outer: DO... inner: DO k=1,5 IF ( a(k)==0 ) CYCLE IF ( a(k)<0 ) EXIT outer... END DO inner END DO outer Exit from loop with EXIT. Begin another iteration with CYCLE. By default EXIT and CYCLE apply to the inner-most loop containing the statement, but can refer to a specific, named loop. INTEGER :: x... IF ( x.lt.1 ) GOTO 10... 10 x = 1
, integer i,j,k,nn integer mm real*8 abc real*8 cde(100),vec(3),xaa(3,3) real*8, allocatable :: dd(:,:), ee(:) integer, allocatable :: kkx(:) : FORTRAN 90 FORTRAN 77.,., FORTRAN 77. allocate( dd(-10:10,3) ) allocate( ee(mm) ). deallocate(dd) ; deallocate(ee),,..,..
C vs FORTRAN F:,. short int x; long int x; int x; char x; char *x; char x[7]; float x; double x; INTEGER*2 X INTEGER*4 X INTEGER*4 X BYTE X or LOGICAL*1 X CHARACTER*n X CHARACTER*7 X REAL X DOUBLE PRECISION X (), C FORTRAN.. a(3,3), a(1,1), a(2,1), a(3,1), a(1,2), a(2,2), a(3,2), a(1,3),a(2,3),a(3,3).,,, b(9)., aaa(100)1,2,3,,,100. C 0,1,2,3.99. FORTRAN bbb(0:99)., FORTRAN, bbb(1:100) bbb(100)., ccc(-100:100).
Slices, section A(:)! the whole array A(3:9)! A(m) to A(n) in steps of 1 A(3:9:1)! as above A(m:n)! A(m) to A(n) A(m:n:k)! A(m) to A(n) in steps of k A(8:3:-1)! A(8) to A(3) in steps of -1 A(8:3)! A(8) to A(3) step 1 => Zero size A(m:)! from A(m) to default UPB A(:n)! from default LWB to A(n) A(::2)! from default LWB to UPB step 2 A(m:m)! 1 element section A(m)! scalar element - not a section aaa(:,3)=bb(:) abc=maxval(aaa(:,2)) abc=minval(abs(aaa(:,:)). (?) FORTRAN 90.
program area implicit none real*8 a,b,c real*8 s,p,aaa read(5,*) a,b,c p=a+b+c s=p/2.0d0 aaa=sqrt(s*(s-a)*(s-b)*(s-c)) write(6,*) aaa stop end program area (5 ), read(5,*) (6 ), write(6,*) a.out < input > output.. implicit none real*8 a,b,c real*8 d read(5,*) a,b,c d=b**2-4.d0*a*c if(d < 0.0d0)then write(6,*) 'NO REAL ROOTS' end if if(d == 0.0d0)then write(6,*) 'two identical roots' write(6,*) -b/(2.d0*a) endif if(d> 0.0d0)then write(6,*) 'two distinct roots' write(6,*) (-b+sqrt(d))/(2.0d0*a) write(6,*) (-b-sqrt(d))/(2.0d0*a) endif stop end
ihlee@cetus0:~/lecture_programming$ a.out Enter nn 10 nn 10 55.00000000000000 aaa ihlee@cetus0:~/lecture_programming$ a.out < input > output., interactive.,.
ihlee@cetus0:~/lecture_programming$ a.out Enter nn 10 nn 10 55.00000000000000 aaa ihlee@cetus0:~/lecture_programming$,... :,,. :,. Subroutine.,.,. Subroutine.
!234567890 program i_sum implicit none integer mm real*8 bbb write(6,*) ' Enter nn' read(5,*) mm write(6,*) 'nn ',mm call int_sum(mm,bbb) write(6,*) bbb,' aaa' end program i_sum subroutine int_sum(nn,aaa) implicit none integer nn real*8 aaa integer i! write(6,*) ' Enter nn'! read(5,*) nn! write(6,*) 'nn ',nn aaa=0.0d0 do i=1,nn aaa=aaa+float(i) end do! write(6,*) aaa,' aaa' end subroutine int_sum ihlee@cetus0:~/lecture_programming$ pgf90 a.f90 sub1.f90 ihlee@cetus0:~/lecture_programming$ a.out Enter nn 100 nn 100 5050.000000000000 aaa ihlee@cetus0:~/lecture_programming$
!234567890 implicit none integer mm real*8 bbb! write(6,*) ' Enter nn' read(5,*) mm write(6,*) 'nn ',mm call int_sum(mm,bbb) write(6,*) bbb,' aaa' call inv_sum(mm,bbb) write(6,*) bbb,' aaa' end subroutine int_sum(nn,aaa) implicit none integer nn real*8 aaa integer i! write(6,*) ' Enter nn'! read(5,*) nn! write(6,*) 'nn ',nn aaa=0.0d0 do i=1,nn aaa=aaa+float(i) end do! write(6,*) aaa,' aaa' end subroutine int_sum ihlee@cetus0:~/lecture_programming$ pgf90 a.f90 sub1.f90 sub2.f90 ihlee@cetus0:~/lecture_programming$ a.out Enter nn 10 nn 10 55.00000000000000 aaa 2.928968253968254 aaa (subroutine, function).
implicit none integer n,nn real*8 a,b real(8) x,f,area,h,s integer k f(x)=x*x! write(6,*) 'a,b,n' read(5,*) a,b,n h=(b-a)/float(n) nn=n-1 s=f(a)+f(b) do k=1,nn s=s+2.0d0*f(a+float(k)*h) enddo area=h*s/2.0d0 write(6,*) a,b,area stop end ihlee@cetus0:~/lecture_programming$ a.out a,b,n 0.d0 1.d0 55 0.0000000000000000 1.000000000000000 0.3333884297520662 FORTRAN STOP ihlee@cetus0:~/lecture_programming$ ihlee@cetus0:~/lecture_programming$ a.out a,b,n 0.d0 1.d0 100 0.0000000000000000 1.000000000000000 0.3333500000000000 FORTRAN STOP ihlee@cetus0:~/lecture_programming$
pgf90 *.f90 pgf90 a.f90 b.f90 c.f90 d.f e.f pgf90 a.f90 b.o c.o d.o e.o makefile
IMPLICIT NONE? IMPLICIT NONE :.,.. REAL*8, INTEGER, CHARACTER*3, LOGICAL, COMPLEX*16. TYPO..,., USE,. IMPLICIT REAL*8 (A-H,O-Z), A-H O-Z., implicit none.... INTEGER, DIMENSION(:), ALLOCATABLE, SAVE, PRIVATE :: ivari_name.,.. INTEGER, ALLOCATABLE :: ivari_name(:) private save,. public private. implicit none private save integer... real*8... real*8, allocatable :: integer, allocatable :: logical... public ::...( ),, ().. Implicit none -.
!234567890 PROGRAM xxx_1 IMPLICIT NONE INTEGER, PARAMETER :: idp=kind(0.d0)! CHARACTER*8 fnnd ; CHARACTER*10 fnnt CHARACTER(8) fnnd ; CHARACTER(10) fnnt LOGICAL ll1,ll2 INTEGER ii,jj,i,nn INTEGER :: kk,npt! REAL*8 aa,pi REAL(8) bb,dx,xx REAL(idp) cc REAL(KIND=idp) dd REAL(KIND=idp) :: ee REAL(8) ff! COMPLEX(idp) zz1,zz2,zz3 COMPLEX(idp) :: zz4! REAL(8), ALLOCATABLE :: xvector(:),xmatrix(:,:) REAL(8), ALLOCATABLE :: yvector(:) REAL(8), ALLOCATABLE :: ymatrix(:,:) REAL(8), ALLOCATABLE :: qvector(:)! REAL(8) :: sma(3) REAL(8) smb(3)! CALL DATE_AND_TIME(DATE=fnnd,TIME=fnnt) WRITE(6,*) ' ',' date ',fnnd,' time ',fnnt,' '! INQUIRE(FILE='input_file',EXIST=ll1)! pi=4.0d0*atan(1.0d0) WRITE(6,*) 'pi ',pi pi=acos(-1.0d0) WRITE(6,*) 'pi ',pi IF(ll1)THEN WRITE(6,*) 'file is present' ELSE WRITE(6,*) 'file is absent' ENDIF WRITE(6,*) idp,' idp' ll1=.true. ll2=.false. WRITE(6,*) ll1,' ll1' WRITE(6,*) ll2,' ll2' ii=0 jj=1 kk=-1 ll1= (1 == 2) ll2= (1 /= 2) IF(ii >= jj) ll1=.true. IF(ii == jj) ll2=.false. WRITE(6,*) ii,' ii' WRITE(6,*) jj,' jj' aa=0.0_idp! note bb=0.0d0 cc=0.0d0 dd=0.0d0 zz1=0.0_idp zz2=cmplx(1.0_idp,2.0_idp)! note zz3=zz1+zz2 zz3=zz1-zz2 zz3=zz1*zz2 zz3=zz1/zz2 zz3=zz1**2 aa=real(zz3) bb=imag(zz3) cc=abs(zz3) zz4=conjg(zz3) WRITE(6,*) aa,' aa' WRITE(6,*) bb,' bb' WRITE(6,*) cc,' cc' WRITE(6,*) dd,' dd' WRITE(6,*) zz1,' zz1' WRITE(6,*) zz2,' zz2' OPEN(1,FILE='output1',FORM='FORMATTED') WRITE(1,'(i8)') ii WRITE(1,'(2i8)') ii,jj WRITE(1,'(i8,2x,i8)') ii,jj WRITE(1,'(i8,2x,i8,3x,a5)') ii,jj,'ii,jj' WRITE(1,*) aa,' aa' WRITE(1,*) bb,' bb' WRITE(1,*) cc,' cc' WRITE(1,*) dd,' dd' WRITE(1,*) zz1,' zz1' WRITE(1,*) zz2,' zz2' CLOSE(1)
OPEN(2,FILE='output2',FORM='FORMATTED') WRITE(2,'(i8)') ii WRITE(2,'(2i8)') ii,jj CLOSE(2) sma=0.0d0 smb=1.0d0 ii=10 jj=10 nn=ii ALLOCATE(xvector(nn)) ; ALLOCATE(yvector(nn)) ALLOCATE(xmatrix(nn,nn)) ALLOCATE(ymatrix(nn,nn)) ALLOCATE(qvector(0:nn)) IF(ALLOCATED(xvector))THEN WRITE(6,*) 'xvector is allocated' WRITE(6,*) 'size ', SIZE(xvector) DO i=1,nn xvector(i)=float(i) ENDDO yvector=-xvector ENDIF WRITE(6,*) 'sum ', SUM(xvector) WRITE(6,*) 'sum ', SUM(yvector) yvector=sqrt(abs(xvector)) IF(ALLOCATED(xmatrix))THEN WRITE(6,*) 'xmatrix is allocated' WRITE(6,*) 'size ', SIZE(xmatrix) ENDIF xvector=0.0d0 yvector=0.0d0 xmatrix=0.0d0 ymatrix=transpose(xmatrix) ymatrix(:,1)=yvector(:) yvector=matmul(xmatrix,xvector) dd=dot_product(xvector,yvector) dd=maxval(xvector) ee=minval(yvector) dd=maxval(abs(xvector)) dd=maxval(abs(xvector-yvector))/(maxval(abs(xvector))+maxval(abs(yvector))+1.0d-8) aa=sum(xvector)/float(size(xvector)) bb=product(yvector) DEALLOCATE(xvector) ; DEALLOCATE(yvector) DEALLOCATE(xmatrix) DEALLOCATE(ymatrix) DEALLOCATE(qvector) INQUIRE(FILE='del',EXIST=ll1) IF(LL1)THEN WRITE(6,*) 'del is present' OPEN(11,FILE='del') CLOSE(31,STATUS='DELETE') ENDIF aa=0.0d0 bb=1.0d0 npt=100000000 dx=(bb-aa)/float(npt-1) cc=0.0d0 DO i=1,npt xx=aa+float(i-1)*dx cc=cc+ff(xx) ENDDO cc=cc*dx write(6,*) cc,' cc' STOP END PROGRAM xxx_1 FUNCTION ff(x) REAL(8) ff REAL(8) x ff=x*x RETURN END FUNCTION ff
ihlee@cetus0:~/lecture_programming$ a.out date 20080109 time 152244.148 pi 3.141592653589793 pi 3.141592653589793 file is absent 8 idp T ll1 F ll2 0 ii 1 jj 0.0000000000000000 aa 0.0000000000000000 bb 0.0000000000000000 cc 0.0000000000000000 dd (0.0000000000000000,0.0000000000000000) zz1 (1.000000000000000,2.000000000000000) zz2 xvector is allocated size 10 sum 55.00000000000000 sum -55.00000000000000 xmatrix is allocated size 100 0.3333333283332683 cc FORTRAN STOP ihlee@cetus0:~/lecture_programming$
20 TOP 10 http://incredible.egloos.com/3216474.,.,.,. Fortran Calls C http://docs.sun.com/app/docs/doc/801-5492/6hvudcqa3?l=ko&a=view C Calls Fortran http://docs.sun.com/app/docs/doc/801-5492/6hvudcqa4?l=ko&a=view
Fortran Calls C void simref_ ( b4, i4, r4, d8 ) /* Simple types, passed by reference, from f90 (f90 calls C)*/ int * b4 ; int * i4 ; float * r4 ; double * d8 ; { *b4 = 1 ; *i4 = 9 ; *r4 = 9.9f ; *d8 = 9.9F ; } PROGRAM SimpleRef! Pass some simple types, by reference, to C (f90 calls C) LOGICAL b4! Default kind is 4-byte INTEGER i4! Default kind is 4-byte REAL r4! Default kind is 4-byte DOUBLE PRECISION d8! This is 8-byte CALL SimRef ( b4, i4, r4, d8 ) WRITE( *, '(L4,I4,F6.1,F6.1)') b4, i4, r4, d8 END PROGRAM SimpleRef ihlee@cetus0:~/lecture_programming$ vi simref.c ihlee@cetus0:~/lecture_programming$ vi simplref.f90 ihlee@cetus0:~/lecture_programming$ cc -c simref.c ihlee@cetus0:~/lecture_programming$ pgf90 simplref.f90 simref.o ihlee@cetus0:~/lecture_programming$ a.out T 9 9.9 9.9 http://docs.sun.com/app/docs/doc/801-5492/6hvudcqa3?l=ko&a=view
,,.,..,.??. : (class) : (object, instance)
-C ( ) -c object. real*8 a(3,2),yone(6) 3*2=6. (row) (column). 1. a(1,1) a(2,1) a(3,1) a(1,2) a(2,2) a(3,2)ij=3*(j-1)+i a(i,j) (yone,, 6.). do i=1,3 do j=1,2 ij=3*(j-1)+i yone(ij)=a(i,j) enddo enddo akvector(natom,3,nk) 1 ( ). do iatom=1,natom kp=3*nk*(iatom-1) do k=1,nk xone(kp+3*(k-1)+1)=akvector(iatom,1,k) xone(kp+3*(k-1)+2)=akvector(iatom,2,k) xone(kp+3*(k-1)+3)=akvector(iatom,3,k) enddo enddo..,, 1..
Subroutine, function( ),. double precision, double precision.. (.) A(3,2). subroutine abc1(a) implicit none real*8 a(3,1) subroutine abc2(a) implicit none real*8 a(3,2).,.,.,. (.) subroutine abc3(a) implicit none real*8 a(1,2) ( 77) subroutine abc4(a,lda) implicit none integer lda real*8 a(lda,1). a(3,3), a(1,1), a(2,1), a(3,1), a(1,2), a(2,2), a(3,2), a(1,3),a(2,3),a(3,3).,,, b(9).
implicit none integer m integer j! read(5,*) m do j=2,m call test_prime(j) end do stop end subroutine test_prime(n) implicit none integer n integer k if( n <=1 ) return if( n == 2) write(6,*) n,' is a prime' if( n == 2) return k=2 30 if (n/k*k == n) goto 70 k=k+1 if(k <= n/2) go to 30 write(6,*) n,' is a prime' return 70 continue write(6,*) n, ' is not a prime' return end ihlee@cetus0:~/lecture_programming$ a.out 20 2 is a prime 3 is a prime 4 is not a prime 5 is a prime 6 is not a prime 7 is a prime 8 is not a prime 9 is not a prime 10 is not a prime 11 is a prime 12 is not a prime 13 is a prime 14 is not a prime 15 is not a prime 16 is not a prime 17 is a prime 18 is not a prime 19 is a prime 20 is not a prime FORTRAN STOP ihlee@cetus0:~/lecture_programming$
!234567890 implicit none integer nm,n,matz,ierr real*8, allocatable :: ar(:,:),ai(:,:),w(:),zr(:,:),zi(:,:),fv1(:),fv2(:),fm1(:,:) real*8 tmp integer i,j,k matz=0 matz=1 nm=3 n=nm open(1,file='input',form='formatted') read(1,*) n close(1) write(6,*) n,' n' allocate(ar(nm,n),ai(nm,n),w(n),zr(nm,n),zi(nm,n),fv1(n),fv2(n),fm1(2,n)) do i=1,n do j=1,i! lower-triangular-form.!------[ set up the hermitian matrix : lower triangular form ar(i,j)=dfloat(i+j)**2 ai(i,j)=1.d0/dfloat(i-j) if(i == j) ai(i,j)=0.0d0 hermitian enddo enddo do i=1,n do j=1,i ar(j,i)=ar(i,j) ai(j,i)=-ai(i,j) enddo enddo call ch(nm,n,ar,ai,w,matz,zr,zi,fv1,fv2,fm1,ierr) write(6,*) 'eigenvalues' do i=1,n write(6,*) w(i) enddo write(6,*) 'eigenvectors' do j=1,n write(6,*) (dcmplx(zr(k,j),zi(k,j)),k=1,n) enddo do i=1,n do j=1,n tmp=0.0d0 do k=1,n tmp=tmp+dcmplx(zr(k,i),-zi(k,i))*dcmplx(zr(k,j),zi(k,j)) enddo write(6,*) i,j,tmp enddo endd deallocate(ar,ai,w,zr,zi,fv1,fv2,fm1)!. stop end 3 3 n eigenvalues -3.295706627065439 0.2569711761519913 59.03873545091345 eigenvectors (0.7426815277007046,0.2134499083118499) (0.3321344158486125,-0.1647469485075425) (-0.5151780835984430,0.000000000000000) (-0.4017627191549581,-0.3764084131040000) (0.7336179945919501,0.2195919189859811) (-0.3323965270105607,0.000000000000000) (0.3152758193792096,-1.9179854283106629E-002) (0.5252646650262870,-1.5040996247447431E-002) (0.7900025892433221,0.000000000000000) 1 1 0.9999999999999998 1 2-2.6099456797157305E-016 1 3-2.2218013019659200E-016 2 1-2.6099456797157305E-016 2 2 0.9999999999999998 2 3-2.0545631168600309E-016 3 1-2.2218013019659200E-016 3 2-2.0545631168600309E-016 3 3 0.9999999999999996
program conjugate_gradient implicit none integer iters, its, n logical :: converged real*8 :: tol, up, alpha, beta real*8, allocatable :: a(:,:), b(:), x(:), r(:), u(:), p(:), xnew(:) read(5,*) n, tol, its allocate( a(n,n), b(n), x(n), r(n), u(n), p(n), xnew(n) ) open(10, file='data') read(10,*) a read(10,*) b close(10) x=1.d0! 1.d0. r=b-matmul(a,x)!.. p=r!. iters=0 do! do loop. iters=iters+1 u=matmul(a,p) up=dot_product(r,r)! (inner product). alpha=up/dot_product(p,u) xnew=x+p*alpha r=r-u*alpha!. beta=dot_product(r,r)/up p=r+p*beta! beta. converged=(maxval(abs(xnew-x))/maxval(abs(x)) < tol)!. x=xnew!. if( converged.or. iters == its) exit! do. enddo write(6,*) iters write(6,*) x deallocate( a,b,x,r,u,p,xnew ) end program conjugate_gradient!.
module common_data implicit none private!.. ( default.) save integer iii1,jjj1 real*, allocatable :: aa1(:,:),bb1(:) logical first_call!.!..true.,.false.. character*2, allocatable :: ccc1(:)... public :: aa1,bb1,iii1!,. (.) end module common_data,.,,,, public. public.,, default public. module common_data2 implicit none private save... public ::... end module common_data2 :. :. Interface.,,. : http://incredible.egloos.com/2950367
( ) CONTAINS (, ). END FUNCTION SUBROUTINE. End function fname,, end subroutine sname... http://www.liv.ac.uk/hpc/htmlhpfcourse/htmlhpfcourseslidesnode142.html
..!234567890 MODULE wf_wr IMPLICIT NONE save public integer nsim real*8 temper real*8, allocatable :: wf(:),wr(:) END MODULE wf_wr!. USE wf_wr IMPLICIT NONE real*8 tol,sol,sol1,sol2 real*8 tmp integer i real*8 bennetteq,zeroin temper=1.d0 ; nsim=10 allocate(wf(nsim),wr(nsim)) wf=1. ; wr=-1. do i=1,nsim call random_number(tmp) wf(i)=wf(i)+0.2*(tmp-0.5) call random_number(tmp) wr(i)=wr(i)+0.2*(tmp-0.5) enddo sol1=maxval(wf) ; sol2=minval(wf) write(6,*) sol1,bennetteq(sol1) write(6,*) sol2,bennetteq(sol2) tol=1.d-16 sol=zeroin(sol1,sol2,bennetteq,tol)! netlib. root finding. : bennetteq. write(6,*) sol,bennetteq(sol) deallocate(wf,wr) stop end
USE, ONLY :, public,. USE USE ONLY.. save, public ( private),. implicit none... 90. implicit none.,,,,,.. -C.. USE stack, local_use_pop => pop pop, stack pop local_pop. USE stack, ONLY: pos, local_pop=> pop stack pos pop., pop local_pop.. makefile., ximage_lbfgs.f90 timpose.f90.. USE aaa USE bbb! bbb public,, (function, subroutine). USE ccc USE ddd, ONLY: qqq,ppp! ddd,, (subroutine, function).
implicit none,, ()..,...,., = {}+{ }. (contain)..,. (, ) USE. 77,, common common. 90 77...,.,.,.,,,, public. public.,, default public.
cmpl=pgf90 segl=pgf90 Ocft= -c -fast cft= -c #LINK=-L/usr/lib/xmm -L/usr/lib/xmm/atlas LINK= -L/usr/lib/sse2 LIBS = -llapack -lf77blas -latlas -lf2c AMD= amd.o predictor_corrector_nose.o energy_force.o tight_binding_lapack_c.o \ mini_driver.o minimization.o lbfgs.o bond_boost.o erf.o maxwell_boltzmann.o \ check_transition.o rmarin.o angle_boost.o dihedral.o dihedral_boost.o amd.x: $(AMD) $(segl) -o amd.x $(AMD) $(LINK) $(LIBS) amd.o:amd.f90 predictor_corrector_nose.o bond_boost.o angle_boost.o dihedral_boost.o $(cmpl) $(Ocft) amd.f90 energy_force.o:energy_force.f90 tight_binding_lapack_c.o $(cmpl) $(Ocft) energy_force.f90 tight_binding_lapack_c.o:tight_binding_lapack_c.f90 $(cmpl) $(Ocft) tight_binding_lapack_c.f90 predictor_corrector_nose.o:predictor_corrector_nose.f90 $(cmpl) $(Ocft) predictor_corrector_nose.f90 mini_driver.o:mini_driver.f90 $(cmpl) $(Ocft) mini_driver.f90 minimization.o:minimization.f90 $(cmpl) $(Ocft) minimization.f90 bond_boost.o:bond_boost.f90 $(cmpl) $(Ocft) bond_boost.f90 angle_boost.o:angle_boost.f90 $(cmpl) $(Ocft) angle_boost.f90 dihedral_boost.o:dihedral_boost.f90 dihedral.o $(cmpl) $(Ocft) dihedral_boost.f90 dihedral.o:dihedral.f90 erf.o:erf.f $(cmpl) $(Ocft) erf.f clean: rm -f *.x *.o *.mod *.M core* touch: touch *.f90 *.i makefile ; chmod 600 *.f90 *.i makefile ; ls -l *.f90 *.i makefile rmo: rm -f *.o *.mod *.M core* lsl: ls -l *.f90 makefile *.i, makefile. make. ls -ltra. ls -Fs make clean make
! module metric implicit none real*8 aa,bb,cc,dd end module metric program abc USE metric implicit none real*8 sol1,sol2,tol real*8 sol,x real*8, external :: func,zeroin integer i aa=0.d0 bb=-2.d0 cc= 1.d0 dd=2.d0 sol1=-1.0d0 sol2= 0.1d0 write(6,*) sol1,func(sol1) write(6,*) sol2,func(sol2) tol=1.d-16 ; sol=zeroin(sol1,sol2,func,tol) write(6,*) sol,func(sol),' sol, f(sol)' sol1=-3.0d0 sol2=-1.9d0 write(6,*) sol1,func(sol1) write(6,*) sol2,func(sol2) tol=1.d-16 ; sol=zeroin(sol1,sol2,func,tol) write(6,*) sol,func(sol),' sol, f(sol)' sol1= 0.9d0 sol2= 1.1d0 write(6,*) sol1,func(sol1) write(6,*) sol2,func(sol2) tol=1.d-16 ; sol=zeroin(sol1,sol2,func,tol) write(6,*) sol,func(sol),' sol, f(sol)' sol1= 1.9d0 sol2= 2.1d0 write(6,*) sol1,func(sol1) write(6,*) sol2,func(sol2) tol=1.d-16 ; sol=zeroin(sol1,sol2,func,tol) write(6,*) sol,func(sol),' sol, f(sol)' stop end program abc real*8 function func(x) USE metric implicit none real*8 x func=(x-dd)*(x-aa)*(x-bb)*(x-cc) return end -1.000000000000000-6.000000000000000 0.1000000000000000 0.3591000000000000 2.7076140240079107E-019 1.0830456096031642E-018 sol, f(sol) -3.000000000000000 60.00000000000000-1.900000000000000-2.148900000000002-2.000000000000000 0.0000000000000000 sol, f(sol) 0.9000000000000000 0.2870999999999999 1.100000000000000-0.3069000000000003 1.000000000000000 0.0000000000000000 sol, f(sol) 1.900000000000000-0.6669000000000005 2.100000000000000 0.9471000000000009 2.000000000000000 0.0000000000000000 sol, f(sol) FORTRAN STOP.,,.,, (.). Cf. 77 common. : http://incredible.egloos.com/3040113
module dihedral_angle! Written by In-Ho Lee, KRISS, October 4 (2003). implicit none private save real*8 diangle,dadi(3),dadj(3),dadk(3),dadl(3) public :: mk_diangle,diangle,dadi,dadj,dadk,dadl contains subroutine mk_diangle(ri,rj,rk,rl) implicit none real*8 ri(3),rj(3),rk(3),rl(3) real*8 rij(3),rkj(3),rkl(3),xnorm_ij_kj,xnorm_kj_kl real*8 rij_kj(3),rkj_kl(3),rtmp(3),xnorm_kj,xnorm_mj,xnorm_nk real*8 rmj(3),rnk(3),xxr,xxi real*8 argu,qsign rij=ri-rj ; rkj=rk-rj ; rkl=rk-rl call cross(rij,rkj,rij_kj) call cross(rkj,rkl,rkj_kl) call cross(rij_kj,rkj_kl,rtmp) xnorm_ij_kj=sqrt(dot_product(rij_kj,rij_kj)) xnorm_kj_kl=sqrt(dot_product(rkj_kl,rkj_kl)) argu=dot_product(rij_kj,rkj_kl)/(xnorm_ij_kj*xnorm_kj_kl) ; argu=min(1.0d0,max(-1.0d0,argu)) qsign=1.0d0 ; if(dot_product(rkj,rtmp) < 0.0d0) qsign=-1.0d0 diangle=qsign*acos(argu) call cross(rij,rkj,rmj) call cross(rkj,rkl,rnk) xnorm_kj=sqrt(dot_product(rkj,rkj)) xnorm_mj=sqrt(dot_product(rmj,rmj)) xnorm_nk=sqrt(dot_product(rnk,rnk)) dadi(:)=rmj(:)*xnorm_kj/xnorm_mj**2 dadl(:)=-rnk(:)*xnorm_kj/xnorm_nk**2 xxr=dot_product(rij,rkj)/xnorm_kj**2-1.0d0 xxi=dot_product(rkl,rkj)/xnorm_kj**2 dadj(:)=xxr*dadi(:)-xxi*dadl(:) xxr=dot_product(rkl,rkj)/xnorm_kj**2-1.0d0 xxi=dot_product(rij,rkj)/xnorm_kj**2 dadk(:)=xxr*dadl(:)-xxi*dadi(:) dadi=dadi*qsign dadj=dadj*qsign dadk=dadk*qsign dadl=dadl*qsign end subroutine mk_diangle subroutine cross(r1,r2,r3) implicit none real*8 r1(3),r2(3),r3(3) r3(1) = r1(2)*r2(3)-r1(3)*r2(2) r3(2) = r1(3)*r2(1)-r1(1)*r2(3) r3(3) = r1(1)*r2(2)-r1(2)*r2(1) end subroutine cross end module dihedral_angle Dihedral angle
!234567890 implicit none real*8, external :: func real*8 rslt,aa,bb integer n n=100 aa=0.0d0 bb=1.0d0 call simpson(func,n,aa,bb,rslt) write(6,*) rslt,' rslt' stop end subroutine simpson(func,n,aa,bb,rslt) implicit none real*8 func integer n real*8 rslt,aa,bb real*8 h,xx integer j logical lodd h=(bb-aa)/float(n) rslt=(func(aa)+func(bb)) lodd=.true. do j=1,n-1 xx=aa+h*float(j) if(lodd)then rslt=rslt+4.0d0*func(xx) else rslt=rslt+2.0d0*func(xx) endif lodd=(.not. lodd) enddo rslt=rslt*h/3.0d0 return end real*8 function func(x) implicit none real*8 x func=x*x return end ---------------------------------------------- 0.3333333333333334 rslt FORTRAN STOP n even number. Composite Simpson's rule : [a,b] step length: h=(b-a)/n x_i=a+i h i=0,1,2,3,...,n-1,n n. x_0=a x_n=b : http://incredible.egloos.com/3021478
, Fornberg ( ) Uniform grid. (nonuniform grid),. (.),. 0 =(k=0) k=1 () k grid. f(x) : http://incredible.egloos.com/3032905
subroutine weights(xi,x,n,m,c)! purpose: computes weights for the mth derivative on arbitrarily! spaced points! Jan 04, 1996 declare all variables! date: Nov 15, 1995 (copied from Fornberg)! written by In-Ho Lee, Beckman Inst., UIUC, April 13 1997.! +real_constant=double option in HP f90! +--------------------------------------------------------------------+! INPUT PARAMETERS:! xi point at which the approximations are to be accurate! x x-coordinates for the grid points, array dimensioned x(0:n)! n the grid points are at x(0),x(1),...x(n) (i.e. n+1 in all)! m highest order of derivative to be approximated!! OUTPUT PARAMETERS! c weights, array dimensioned c(0:n,0:n,0:m)! on return, the element c(j,k,i) contains the weight to be! applied at x(k) when the i:th derivative is approximated! by a stencil extending over x(0),x(1),...,x(j).! +--------------------------------------------------------------------+! implicit real*8 (a-h,o-z) dimension x(0:n),c(0:n,0:n,0:m) c(0,0,0)= 1.0d0 c1 = 1.0d0 c4 = x(0)-xi do 40 j=1,n mn = min(j,m) c2 = 1.0d0 c5 = c4 c4 = x(j)-xi do 20 k=0,j-1 c3=x(j)-x(k) c2=c2*c3 if (j.le.m) c(k,j-1,j)=0.0d0 c(k,j,0) =c4*c(k,j-1,0)/c3 do 10 i=1,mn 10 c(k,j,i)=(c4*c(k,j-1,i)-i*c(k,j-1,i-1))/c3 20 continue c(j,j,0)=-c1*c5*c(j-1,j-1,0)/c2 do 30 i=1,mn 30 c(j,j,i)=c1*(i*c(j-1,j-1,i-1)-c5*c(j-1,j-1,i))/c2 40 c1=c2 return end, Fornberg ( ) Fornberg Finite difference., high-order.. Finite-difference. 10.,. (, : ).,, m.,,... Finite difference: http://en.wikipedia.org/wiki/finite_difference http://mathworld.wolfram.com/finitedifference.html http://en.wikipedia.org/wiki/numerical_differentiation
!234567890 module rk_df implicit none public save integer neq,info(15),idid,lrw,liw real*8 t, tout integer, allocatable :: iwork(:),ipar(:) real*8, allocatable :: y(:), rwork(:), rtol(:), atol(:), rpar(:) contains subroutine initial implicit none t=0.0d0 tout=0.1d0 neq=3 allocate(y(neq)) y(1)=0.d0 y(2)=1.d0 y(3)=1.d0 info(1)=0 lrw= 33+7*neq ; liw= 34 allocate(iwork(liw)) ; allocate(rwork(lrw)) allocate(atol(neq),rtol(neq)) rtol=1.d-10 ; atol=1.d-10 allocate(rpar(neq)) ; allocate(ipar(neq)) end subroutine initial subroutine final implicit none deallocate(y) ; deallocate(rwork) ; deallocate(iwork) deallocate(atol,rtol) ; deallocate(rpar) ; deallocate(ipar) end subroutine final! DF(X,U,UPRIME,RPAR,IPAR) x, u.! uprime(). subroutine df(x,u,uprime,rpar_z,ipar_z) implicit none real*8 x,u(1) real*8 rpar_z(1) integer ipar_z(1) real*8 uprime(1) uprime(1)=u(2)*u(3) uprime(2)=-u(1)*u(3) uprime(3)=-0.522d0*u(1)*u(2) end subroutine df end module rk_df! program rk_test USE rk_df, ONLY : neq, t, y, tout, info, rtol, atol, idid, rwork, lrw, iwork, liw, rpar, ipar USE rk_df, ONLY : df,initial,final implicit none call initial call dderkf(df, neq, t, y, tout, info, rtol, atol, idid, rwork, lrw, iwork, liw, rpar, ipar)! write(6,*) t,' t' write(6,*) tout,' tout' write(6,*) y,' y' write(6,'(i5,1x,a4)') idid,'idid' write(6,'(15i5,1x,a4)') info,'info' tout=tout+0.1d0 call dderkf(df, neq, t, y, tout, info, rtol, atol, idid, rwork, lrw, iwork, liw, rpar, ipar)! write(6,*) t,' t' write(6,*) tout,' tout' write(6,*) y,' y' write(6,'(i5,1x,a4)') idid,'idid' write(6,'(15i5,1x,a4)') info,'info' tout=tout+0.1d0 call dderkf(df, neq, t, y, tout, info, rtol, atol, idid, rwork, lrw, iwork, liw, rpar, ipar)! write(6,*) t,' t' write(6,*) tout,' tout' write(6,*) y,' y' write(6,'(i5,1x,a4)') idid,'idid' write(6,'(15i5,1x,a4)') info,'info' call final end program rk_test allocate(vector(-89:100)). -89, -88, -87,...98, 99, 100., allocate(vector(100)) 100 vector., 1 100. allocate(vector(1:100)). 1. : http://incredible.egloos.com/3143433
implicit none integer j,i integer narr,mref real*8, allocatable :: xarr(:,:) real*8, allocatable :: yarr(:,:) logical lexist narr=10 mref=3 allocate(xarr(narr,mref)) allocate(yarr(narr,mref)) inquire(file='fort.2',exist=lexist) if(lexist)then open(2,file='fort.2',access='direct',recl=8*narr,form='unformatted') do j=1,mref read(2,rec=j) yarr(:,j) enddo close(2) endif do j=1,mref do i=1,narr xarr(i,j)=float(i)**2/float(j) enddo enddo write(6,*) maxval(abs(xarr-yarr)) open(2,file='fort.2',access='direct',recl=8*narr,form='unformatted') do j=1,mref write(2,rec=j) xarr(:,j) enddo close(2) deallocate(xarr) deallocate(yarr) stop end implicit none integer j,i integer narr,mref real*8, allocatable :: xarr(:,:) real*8, allocatable :: yarr(:,:) logical lexist narr=10 mref=3 allocate(xarr(narr,mref)) allocate(yarr(narr,mref)) inquire(file='fort.12',exist=lexist) if(lexist)then open(12,file='fort.12',form='unformatted') read(12) yarr(:,:) close(12) endif do j=1,mref do i=1,narr xarr(i,j)=float(i)**2/float(j) enddo enddo write(6,*) maxval(abs(xarr-yarr)) open(12,file='fort.12',form='unformatted') write(12) xarr(:,:) close(12) deallocate(xarr) deallocate(yarr) stop end : http://incredible.egloos.com/3044815
READ(*,100) i,j WRITE(*,100) i,j READ(*,FMT=200) x,y WRITE(*,200) x,y 100 FORMAT(2I) 200 FORMAT(2F10.6) READ(*,'(2I)') i,j WRITE(*,'(2F12.6)') x,y WRITE(6,`(2F12.6)`) 12.6, -131.45678 INTEGER :: i, j REAL, DIMENSION(10) :: a REAL, DIMENSION(10,10) :: b READ(*,*) (a(j),j=1,10) WRITE(*,*) (a(j), j=10,1,-1) WRITE(*,*) ((b(i,j),i=1,10), j=1,10). formatted, unformatted. Unformatted. /... form= formatted serial, direct. Direct CD, 3 7. access= direct,. recl. Record length.
open(11,file= abcd.txt, form= unformatted, access= direct, recl=8*nnsszz) Unformatted... write(11,rec=k) (aa(i),i=1,nnsszz) Unit number... read(11,rec=kk) (bb(j),j=1,nnsszz) /(record length). Unformatted read/write... cf. netcdf, HDF hdf.ncsa.uiuc.edu/products/hd5/index.html
Inquiry list EXIST=lex!.true. or.false. OPENED=lod!.true. or.false. NUMBER=unum! unit number NAME=fnm! filename ACCESS=acc! 'DIRECT' or 'SEQUENTIAL' SEQUENTIAL=seq! 'YES or 'NO' DIRECT=dir! 'YES' or 'NO' FORMATTED=fmt! 'YES' or 'NO' UNFORMATTED=unfmt! 'YES' or 'NO' FORM=frm! 'FORMATTED' or 'UNFORMATTED' NEXTREC=recn! number of next record RECL=recl! record length inquire(file='fort.40',exist=lexist) if(lexist)then else. endif allocated(abc) abc..true.,.false.. if(allocated(abc)) deallocate(abc) if(.not. allocated(abc)) allocate(abc(0:100)).
implicit none real*8, allocatable :: xdata(:),ydata(:) character*15, file_a,file_b integer j,jtemp,ndata real*8 rr,prob,zz write(6,*) 'Enter: ndata' read(5,*) ndata write(6,*) 'Enter: file_a' read(5,*) file_a write(6,*) file_a write(6,*) 'Enter: file_b' read(5,*) file_b write(6,*) file_b allocate(xdata(ndata),ydata(ndata)) open(11,file=file_a,form='formatted') do j=1,ndata read(11,*) jtemp,xdata(j) enddo close(11) open(21,file=file_b,form='formatted') do j=1,ndata read(21,*) jtemp,ydata(j) enddo close(21) call pearsn(xdata,ydata,ndata,rr,prob,zz) write(6,*) rr,rr**2, ' r, r^2' write(6,*) prob, ' prob' write(6,*) zz, ' z' deallocate(xdata,ydata) STOP END Linear correlation -1.1 : http://incredible.egloos.com/1992323
implicit none integer j, k,i real*8 tmp real cpu_diff real cpu_time1 real cpu_time2 real ( kind = 8 ) real_diff real ( kind = 8 ) real_time1 real ( kind = 8 ) real_time2 call cpu_time (cpu_time1 )! call real_time ( real_time1 ) tmp=0.0d0 do i=1,1000 do j=1,1000 do k=1,100 tmp=tmp+sin(float(j*i))*cos(float(i*j)) enddo enddo enddo ihlee@cetus0:~/lecture_programming$ a.out 16.00000 FORTRAN STOP ihlee@cetus0:~/lecture_programming$ ihlee@cetus0:~/lecture_programming$ time a.out 16.00000 FORTRAN STOP real 0m16.031s user 0m16.028s sys 0m0.001s ihlee@cetus0:~/lecture_programming$ call cpu_time ( cpu_time2 )! call real_time ( real_time2 ) cpu_diff = cpu_time2 - cpu_time1! real_diff = real_time2 - real_time1 write(6,*) cpu_diff! write(6,*) real_diff stop end : http://incredible.egloos.com/3047318
subroutine timestamp ( )!*****************************************************************************80!!! TIMESTAMP prints the current YMDHMS date as a time stamp.!!example:!! May 31 2001 9:45:54.872 AM!! Modified:!! 31 May 2001!! Author:!! John Burkardt!! Parameters:!! None! implicit none character ( len = 8 ) ampm integer d character ( len = 8 ) date integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s character ( len = 10 ) time integer values(8) integer y character ( len = 5 ) zone call date_and_time ( date, time, zone, values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0.and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0.and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if!234567890 implicit none call timestamp() stop end ----------------------------------------- March 12 2007 1:58:40.868 PM FORTRAN STOP write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end : http://incredible.egloos.com/3047077
subroutine lbfgs_mini(object,posi,grad,nn)! Written by In-Ho Lee, KRISS, March 3 (2003). implicit none integer nn real*8 object,posi(nn),grad(nn) integer n_local,m,msave real*8, allocatable :: diag(:),w(:) real*8 eps,xtol,gtol,t1,t2,stpmin,stpmax integer iprint(2),iflag,icall,mp,lp,nwork integer jnatom,j logical diagco real*8, allocatable :: xyz(:,:),fxyz(:,:)! The driver for LBFGS must always declare LB2 as EXTERNAL external lb2 common /lb3/mp,lp,gtol,stpmin,stpmax n_local=nn m=5 ; msave=7 ; nwork=n_local*(2*msave+1)+2*msave allocate(diag(n_local),w(nwork)) iprint(1)= 1 ; iprint(2)= 0!! We do not wish to provide the diagonal matrices Hk0, and! therefore set DIAGCO to FALSE. diagco=.false. ; eps= 1.0d-4 ; xtol= 1.0d-16 ; icall=0 ; iflag=0! jnatom=nn/3 allocate(xyz(jnatom,3),fxyz(jnatom,3))! 20 continue do j=1,jnatom xyz(j,1)=posi(3*(j-1)+1) xyz(j,2)=posi(3*(j-1)+2) xyz(j,3)=posi(3*(j-1)+3) enddo call energy_force(xyz,fxyz,object,jnatom) do j=1,jnatom grad(3*(j-1)+1)=-fxyz(j,1) grad(3*(j-1)+2)=-fxyz(j,2) grad(3*(j-1)+3)=-fxyz(j,3) enddo! call lbfgs(n_local,m,posi,object,grad,diagco,diag,iprint,eps,xtol,w,iflag) if(iflag <= 0) go to 50 icall=icall + 1! We allow at most 2000 evaluations of Function and Gradient if(icall > 2000) go to 50 go to 20 50 continue deallocate(diag,w) deallocate(xyz,fxyz) return end ( ) 2 1 lbfgs 1. 2. Gradient, function : : http://incredible.egloos.com/3042973
!234567890 implicit none integer nm,n,matz,ierr real*8, allocatable :: ar(:,:),ai(:,:),w(:),zr(:,:),zi(:,:),fv1(:),fv2(:),fm1(:,:) real*8 tmp integer i,j,k matz=0 matz=1 nm=3 n=nm open(1,file='input',form='formatted') read(1,*) n close(1) write(6,*) n,' n' allocate(ar(nm,n),ai(nm,n),w(n),zr(nm,n),zi(nm,n),fv1(n),fv2(n),fm1(2,n)) do i=1,n do j=1,i! lower-triangular-form.!------[ set up the hermitian matrix : lower triangular form ar(i,j)=dfloat(i+j)**2 ai(i,j)=1.d0/dfloat(i-j)!------] if(i == j) ai(i,j)=0.0d0 hermitian enddo enddo do i=1,n do j=1,i ar(j,i)=ar(i,j) ai(j,i)=-ai(i,j) enddo enddo call ch(nm,n,ar,ai,w,matz,zr,zi,fv1,fv2,fm1,ierr) write(6,*) 'eigenvalues' do i=1,n write(6,*) w(i) enddo write(6,*) 'eigenvectors' do j=1,n write(6,*) (dcmplx(zr(k,j),zi(k,j)),k=1,n) enddo do i=1,n do j=1,n tmp=0.0d0 do k=1,n tmp=tmp+dcmplx(zr(k,i),-zi(k,i))*dcmplx(zr(k,j),zi(k,j)) enddo write(6,*) i,j,tmp enddo enddo deallocate(ar,ai,w,zr,zi,fv1,fv2,fm1)!. stop end Matrix diagonalization( )
!234567890 implicit none integer i,j,k character*32 filename character*4 fn integer nsize nsize=4 do j=0,3 call numeral (j,fn,nsize) write(6,*) fn filename=' del'//fn write(6,*) filename! call system("cp admd.seq"//filename) end do stop end ihlee@cetus0:~/lecture_programming$ a.out 0000 del0000 0001 del0001 0002 del0002 0003 del0003 FORTRAN STOP ihlee@cetus0:~/lecture_programming$ : http://incredible.egloos.com/3474225
!234567890 implicit none integer i,j,k character*32 subcommand character*32 subcommand2 character*4 fn character*9 fm integer nsize nsize=4 do j=0,3 call numeral (j,fn,nsize) write(6,*) fn fm=' admd'//fn subcommand=fm//'.pdb' write(6,*) subcommand! call system("cp admd.pdb"//subcommand) subcommand2=fm//'.dssp' subcommand2=trim(subcommand)//fm//'.dssp' write(6,*) subcommand2! call system("~/dssp/dsspcmbi "//subcommand2) end do stop end : http://incredible.egloos.com/3474225
MODULE cartesian TYPE point REAL :: x, y END TYPE point INTERFACE OPERATOR (.DIST. ) MODULE PROCEDURE dist END INTERFACE CONTAINS REAL FUNCTION dist( a, b ) TYPE(point), INTENT(IN) :: a, b dist = SQRT( (a%x-b%x)**2 + (a%y-b%y)**2 ) END FUNCTION dist END MODULE cartesian ihlee@cetus0:~/lecture_programming$ a.out 1.414214 FORTRAN STOP ihlee@cetus0:~/lecture_programming$ program over2 USE cartesian IMPLICIT NONE TYPE(point) :: a, b REAL distance a%x=0.0d0 a%y=0.0d0 b%x=1.0d0 b%y=1.0d0 distance = a.dist. b write(6,*) distance stop end program over2
MODULE cartesian TYPE point REAL*8 :: x, y, z END TYPE point INTERFACE OPERATOR (.DIST. ) MODULE PROCEDURE dist END INTERFACE CONTAINS REAL*8 FUNCTION dist( a, b ) TYPE(point), INTENT(IN) :: a, b dist = SQRT( (a%x-b%x)**2 + (a%y-b%y)**2 + (a%z-b%z)**2 ) END FUNCTION dist END MODULE cartesian program over2 USE cartesian IMPLICIT NONE TYPE(point) :: a, b REAL*8 distance a%x=0.0d0 a%y=0.0d0 a%z=0.0d0 b%x=1.0d0 b%y=1.0d0 b%z=1.0d0 ihlee@cetus0:~/lecture_programming$ a.out 1.732050807568877 FORTRAN STOP ihlee@cetus0:~/lecture_programming$ TYPE COORDS_3D REAL :: x, y, z END TYPE COORDS_3D TYPE(COORDS_3D) :: pt1, pt2 distance = a.dist. b write(6,*) distance stop end program over2
MODULE cartesian TYPE point REAL*8 :: x, y END TYPE point INTERFACE ASSIGNMENT( = ) MODULE PROCEDURE max_point END INTERFACE CONTAINS SUBROUTINE max_point( a, pt ) REAL*8, INTENT(OUT) :: a TYPE(point), INTENT(IN) :: pt a = MAX( pt%x, pt%y ) END SUBROUTINE max_point END MODULE cartesian ihlee@cetus0:~/lecture_programming$ a.out 4.200000000000000 FORTRAN STOP ihlee@cetus0:~/lecture_programming$ program over4 USE cartesian IMPLICIT NONE TYPE(point) :: a = point(1.7d0, 4.2d0) REAL*8 :: coord coord = a write(6,*) coord stop end program over4
PROGRAM test INTERFACE REAL*8 FUNCTION ratio(x, y) REAL*8 ::x, y END FUNCTION ratio END INTERFACE ihlee@cetus0:~/lecture_programming$ a.out The ratio is 0.1200000000000000 ihlee@cetus0:~/lecture_programming$ INTEGER :: i=3, j=25 PRINT *,'The ratio is ',ratio(dfloat(i),dfloat(j)) END PROGRAM test REAL*8 FUNCTION ratio(x, y) REAL*8 :: x, y ratio=x/y END FUNCTION ratio
PROGRAM test!interface! REAL*8 FUNCTION ratio(x, y)! REAL*8 ::x, y! END FUNCTION ratio!end INTERFACE ihlee@cetus0:~/lecture_programming$ a.out The ratio is 0.1200000000000000 ihlee@cetus0:~/lecture_programming$ INTEGER :: i=3, j=25 real*8 ratio PRINT *,'The ratio is ',ratio(dfloat(i),dfloat(j)) END PROGRAM test REAL*8 FUNCTION ratio(x, y) REAL*8 :: x, y ratio=x/y END FUNCTION ratio Procedure interface.
SUBROUTINE invert(a, inverse, count) REAL, INTENT(IN) :: a REAL, INTENT(OUT) :: inverse INTEGER, INTENT(INOUT) :: count inverse = 1/a count = count+1 END SUBROUTINE invert PROGRAM test INTERFACE REAL FUNCTION func( x ) REAL, INTENT(IN) ::x END FUNCTION func END INTERFACE... CALL sub1( a, b, func(2) )... END PROGRAM test REAL FUNCTION func( x )! external REAL, INTENT(IN) :: x func = 1/x END FUNCTION func
MODULE strings INTERFACE OPERATOR ( / ) MODULE PROCEDURE num END INTERFACE CONTAINS INTEGER FUNCTION num( s, c ) CHARACTER(len=*), INTENT(IN) :: s CHARACTER, INTENT(IN) :: c num = 0 DO i=1,len( s ) IF( s(i:i)==c ) num=num+1 END DO END FUNCTION num END MODULE strings program ov1 USE strings implicit none character*32, ch1,ch2,ch3 integer ii,jj,kk ihlee@cetus0:~/lecture_programming$ a.out 3 1 1 FORTRAN STOP ihlee@cetus0:~/lecture_programming$ Operator overloading ch1='hello world' ch2='my name is in-ho lee' ch3='who are you' ii = ch1/'l' jj = ch2/'o' kk = ch3/'a' write(6,*) ii write(6,*) jj write(6,*) kk stop end program ov1 http://www.qmw.ac.uk/~cgaa260/building/intr_f90/prg_unts/opoload.htm
MODULE SWAPPER INTERFACE SWAP MODULE PROCEDURE SWAP_R, SWAP_I, SWAP_C END INTERFACE CONTAINS SUBROUTINE SWAP_R(A, B) IMPLICIT NONE REAL, INTENT (INOUT) :: A, B REAL :: TEMP TEMP = A ; A = B ; B = TEMP END SUBROUTINE SWAP_R ihlee@cetus0:~/lecture_programming$ a.out 1 2 100 200 7.100000 10.90000 11.10000 17.00000 ab1" 2 1 200 100 10.90000 7.100000 17.00000 11.10000 ba"1 ihlee@cetus0:~/lecture_programming$ SUBROUTINE SWAP_I(A, B) IMPLICIT NONE INTEGER, INTENT (INOUT) :: A, B INTEGER :: TEMP TEMP = A ; A = B ; B = TEMP END SUBROUTINE SWAP_I SUBROUTINE SWAP_C(A, B) IMPLICIT NONE CHARACTER, INTENT (INOUT) :: A, B CHARACTER :: TEMP TEMP = A ; A = B ; B = TEMP END SUBROUTINE SWAP_C END MODULE SWAPPER PROGRAM SWAP_MAIN USE SWAPPER IMPLICIT NONE INTEGER :: I, J, K, L REAL :: A, B, X, Y CHARACTER :: C, D, E, F I = 1 ; J = 2 ; K = 100 ; L = 200 A = 7.1 ; B = 10.9 ; X = 11.1; Y = 17.0 C = 'a' ; d = 'b' ; E = '1' ; F = '"' WRITE (*,*) I, J, K, L, A, B, X, Y, C, D, E, F CALL SWAP (I, J) ; CALL SWAP (K, L) CALL SWAP (A, B) ; CALL SWAP (X, Y) CALL SWAP (C, D) ; CALL SWAP (E, F) WRITE (*,*) I, J, K, L, A, B, X, Y, C, D, E, F END http://www-teaching.physics.ox.ac.uk/unix+prog/hargrove/tutorial_90/13_extra.html
Derived Type TYPE COORDS_3D REAL :: x, y, z END TYPE COORDS_3D TYPE(COORDS_3D) :: pt1, pt2 TYPE SPHERE TYPE (COORDS_3D) :: centre REAL :: radius END TYPE SPHERE TYPE (SPHERE) :: bubble, ball pt1%x = 1.0 bubble%radius = 3.0 bubble%centre%x = 1.0 pt1 = COORDS_3D(1.,2.,3.) bubble%centre = COORDS_3D(1.,2.,3.) bubble = SPHERE(bubble%centre,10.) bubble = SPHERE(COORDS_3D(1.,2.,3.),10.)
Multiple Precision Computation. 2**200 = 1606938044258990275541962092341162602522202993782792835301376 (prime number)? 5468317884572019103692012212053793153845065543480825746529998049913561 The same example using derived types:... USE FMZM TYPE ( FM ) A,B,C,D... D = SQRT( A**2 + B**2 + C**2 )... http://incredible.egloos.com/3160240
References http://incredible.egloos.com/2950367 pdf.. 90. http://www.fortran.com/fortran/ Pointer to everything Fortran http://www.arc.unm.edu/~jabed/list.html Pointer to a list of tutorials http://dynaweb.sdsc.edu:8080/library/t90_c90_y90_cray_fp?dweb_collection=dweb_collections Cray's Manual http://www.software.ibm.com/ad/fortran/xlfortran/optim.html Short optimization guide from IBM http://www.digital.com/fortran/dvf.html DEC's Digital Visual Fortran Page ftp://mycroft.plk.af.mil/pub/fortran_90/tutorial A tutorial by Zane Dodson http://www.nova.edu/ocean/psplot.html Postscript plotting library http://meteora.ucsd.edu/~pierce/fxdr_home_page.html Subroutines to do unformatted I/O across platforms, provided by David Pierce at UCSD http://www.nsc.liu.se/~boein/f77to90/a5.html A good reference for intrinsic functions http://www.nsc.liu.se/~boein/f77to90/ Fortran 90 for the Fortran 77 Programmer http://www.tat.physik.uni-tuebingen.de/~kley/lehre/ftn77/tutorial/index.html http://www.sdsc.edu/~tkaiser/f90.html
(?)
Very-high-level language (shell) (very-high-level language). (.). (!). (; ; object-oriented)..... No compile Monty Python's Flying Circus. No type declarations No memory allocation control High level data structure : (Guido van Rossum) http://en.wikibooks.org/wiki/python_programming http://incredible.egloos.com/2950369 http://en.wikipedia.org/wiki/python_(programming_language) http://pentangle.net/python/handbook/handbook.html http://www.pasteur.fr/recherche/unites/sis/formation/python/
. (dynamic typing),., (),,. (glue language)..,. (dynamic typing). (.).,,. ("Battery Included"),.. (MS windows,,,.,.)! : http://www.awaretek.com/tutorials.html#begin http://ko.wikipedia.org/wiki/%ed%8c%8c%ec%9d%b4%ec%8d%ac
Python Control-d exit(). ihlee@cetus0:~/lecture_programming$ python Python 2.3.5 (#2, Oct 16 2006, 19:59:54) [GCC 3.3.5 (Debian 1:3.3.5-13)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> print "Hello world" ihlee@cetus0:~/lecture_programming$ python Python 2.3.5 (#2, Oct 16 2006, 19:59:54) [GCC 3.3.5 (Debian 1:3.3.5-13)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> a=1. ; b=2.d0 File "<stdin>", line 1 a=1. ; b=2.d0 ^ SyntaxError: invalid syntax >>> a=1 ; b=2. >>> c=a+b >>> print c 3.0 >>> bc l. quit. ihlee@cetus0:~$ bc -l bc 1.06 Copyright 1991-1994, 1997, 1998, 2000 Free Software Foundation, Inc. This is free software with ABSOLUTELY NO WARRANTY. For details type `warranty'. 12.*12 144 quit ihlee@cetus0:~$
* (indentation) python.! *# comment. ``. * \.) "". * (,,,. *. #!/usr/bin/env python. * "def" ":".. if for. *for i in range(100): i 0 0,1,2,3,4,...99 for loop 100. range range(1,101) 1 100 i. range(10) ( ). [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] *FORTRAN 90.. * return.,... return. return a,b,c,d. a, b 1, c 2 ( ) d.. abcd() aa,bb,cc,dd=abcd().
#!/usr/bin/env python #from Numeric import * #from string import atof,atoi,split #from posix import getpid,rename,remove,mkdir,rmdir,system #from shutil import copyfile #import os,sys,time,random version ="1.0.0" counter=0 acc = 0 for n in xrange(1,10000000000): num1 = str(n).count('5') if num1 > 0 : acc = acc+num1 if acc == n: counter=counter+1 print "counter = ",counter,n if counter == 2 : break,. chmod 755 abc.py
system( )... ` `;? CALL SYSTEM('application.py'). call system. ( ).... : / == call system('python_driver.py') ==. : / ==
sad.x 2000.,.. for i in range(2000): # 2000 gen_sad_input() # sad.i system('sad.x >TMPFILE') # sad.x copyfile('dim.dat', ascnum('dimer.',i)) copyfile('mid.xyz', ascnum_xyz('midpt.',i)) # system('rm dim.dat') # system('rm TMPFILE') os.path.exists >>> from os import * >>> system('ls') abc.py HFB PBS_script_examples tr.py action_protein_draco JOB_ID_PRINTED swkim_mpi VASP_example admd_cluster matplot_test test v_tubes dssp onion test1 ykm 0 >>> system(), system('qsub work.0003/pbs_script_file') pbs. j=3 ; f_name=ascnum('work.',j)+'/pbs_script_file' ; system('qsub '+f_name+'') import glob,time import os.path file_list=glob.glob('*') for f_name in file_list: if os.path.isfile(f_name): print f_name,' is a regular file' print os.path.getsize(f_name) print os.path.getatime(f_name) print time.ctime(os.path.getatime(f_name)) elif os.path.isdir(f_name): print f_name,' is a directory' elif os.path.islink(f_name): print f_name,' is a symbolic link'
#!/usr/bin/env python faren_old = float(raw_input("enter yesterday's temperature ")) faren_new = float(raw_input("enter today's temperature ")) if faren_new > 90 and faren_old < 90: print "First day of a heatwave" elif faren_new > 90: print "Subsequent day of a heatwave" elif faren_old < 90: print "Not the hot season" else: print "Bye, bye heatwave" print "Program completed" >>> raw_input() 3 '3' >>> input() 4 4 ihlee@cetus0:~/lecture_programming$ py1.py Enter yesterday's temperature 99 Enter today's temperature 98. Subsequent day of a heatwave Program completed ihlee@cetus0:~/lecture_programming$
phi_trj=[[0.]*(iidd) for j in xrange(np+1)] psi_trj=[[0.]*(iidd) for j in xrange(np+1)] for ifile in range(0,np+1): tmp_file=asc2num('rama',ifile)+'.txt' fp_rama=open(tmp_file,'r') jth=0 while 1: line=fp_rama.readline() if not line : break else: token=line.split() phi_trj[ifile][jth]=float(token[0]) psi_trj[ifile][jth]=float(token[1]) jth=jth+1 fp_rama.close() http://incredible.egloos.com/3588647
def wwait(file_name): while 1 : time.sleep(60.) if os.path.exists(file_name) : remove(file_name) break return # 60 # def mkdir_series(n_dir): for j in range(n_dir): char=ascnum('test_dir.',j) # test_dir.0001 n_dir mkdir(char) return touch abcd
#!/usr/bin/env python from Numeric import * from string import atof,atoi,split from posix import getpid,rename,remove,mkdir,rmdir,system from shutil import copyfile import os,sys,time,random version ="1.0.0" if name == " main ": # # *.bin # convert.x #. convert.x "convert.x test.bin test". # # yeskim98, March/26/2003 # system('ls *.bin >namelist') # f=open('namelist','r') for line in f: name=line.split()[0] name1=name[:-4] # 4 system('./convert.x '+name+' '+name1+'') f.close()
#!/usr/bin/env python from Numeric import * # from string import atof,atoi,split from posix import getpid,rename,remove,mkdir,rmdir,system from shutil import copyfile import os,sys,time,random version ="1.0.0" def prtjobid(): # written by In-Ho Lee, KRISS, Jan. 12, 2002. f=open('job_id_printed','w') f.write('%15d \n' % os.getpid()) # st=time.asctime() # f.write('%23s \n' % st) f.close() return def num2char(i): # ic=str(i) ic='0' if i == 1 : ic='1' if i == 2 : ic='2' if i == 3 : ic='3' if i == 4 : ic='4' if i == 5 : ic='5' if i == 6 : ic='6' if i == 7 : ic='7' if i == 8 : ic='8' if i == 9 : ic='9' return ic def ascnum(ichar,number): # written by In-Ho Lee, KRISS, Jan. 12, 2002. i=int(float(number)/1000.) j=int((number-i*1000.)/100.) k=int((number-j*100.-i*1000.)/10.) l=int(number-k*10.-j*100.-i*1000.) ic=num2char(i) jc=num2char(j) kc=num2char(k) lc=num2char(l) fname=ichar+ic+jc+kc+lc # abc.0001 return fname ihlee@cetus0:~/lecture_programming$ num.py abc0000.xyz abc0001.xyz abc0002.xyz abc0003.xyz abc0004.xyz abc0005.xyz abc0006.xyz abc0007.xyz abc0008.xyz abc0009.xyz ihlee@cetus0:~/lecture_programming$ >>> a=1 >>> b=str(a).zfill(4) >>> print b 0001 >>> c=str(a).zfill(5) >>> print c 00001 >>> c=str(c) >>> print c 00001 >>> >>> from string import * >>> a=1 >>> b=str(a) >>> print b 1 >>> a=1.e0 >>> c=str(222) >>> a=complex(a) >>> print c >>> print a 222 (1+0j) >>> >>> def ascnum_xyz(ichar,number): # written by In-Ho Lee, KRISS, Jan. 12, 2002. fname=ascnum(ichar,number)+'.xyz' # abc.0001.xyz return fname prtjobid() ichar='abc' for i in range(10): print ascnum_xyz(ichar,i) >>> print atof(c) 222.0 >>> print atoi(c) 222 >>> print float(c) 222.0 >>> print int(c) 222
Prime number? #!/usr/bin/env python for i in range(2,10): for x in range(2,i): if i % x == 0 : print i, 'equals', x, '*', i/x break else: print i, 'is a prime number' ihlee@cetus0:~/lecture_programming$ prime.py 3 is a prime number 4 equals 2 * 2 5 is a prime number 5 is a prime number 5 is a prime number 6 equals 2 * 3 7 is a prime number 7 is a prime number 7 is a prime number 7 is a prime number 7 is a prime number 8 equals 2 * 4 9 is a prime number 9 equals 3 * 3 ihlee@cetus0:~/lecture_programming$
a.lt. b 77 a.le. b 77 a.gt. b a.ge. b a.eq. b a /= b 90 a.ne. b a < b.and. b < c a < b a <= b a > b a >= b a == b a!= b a < b < c a is less than b a is less than or equal to b a is greater than b a is greater than or equal to b a is equal to b a is not equal to b a is less than b, which is less than c 90, python =. >=, <=,!=(/=).
if x == 10 and y > z: print "Some statements which only get executed if" print "x is equal to 10 AND y is greater than z." if x == 10 or y > z: print "Some statements which get executed if either print "x is equal to 10 OR y is greater than z" if [condition]: [Some statements executed only if [condition] is true] else: [Some statements executed only if [condition] is false] if 0 <= x <= 10: print "That is between zero and ten inclusive" elif 10 < x < 20: print "That is between ten and twenty" else: print "That is outside the range zero to twenty" FORTRAN 90 endif. elif.. : 1 : 0
Reserved Words and assert break class continue def del elif else except exec finally for from global if import in is lambda not or pass print raise return try while Data Float Int Numeric Oxphys array close float int input open range type write zeros acos asin atan cos e exp fabs floor log log10 pi sin sqrt tan reserved word., Numeric.py,..
FORTRAN 90 :.and. ======= and.or. ======= or.gt. (>) ======= >.ge. (>=) ======= >= Shortcut Operators.le. (<=) ======= <=.eq. (==) ======== ==.ne. (/=) ========!= a < b < c a is less than b, which is less than c 2 2**2. %. not. not logical_variable. sqrt() exp() log() log10() sin() cos() tan() asin() acos() atan() floor() fabs() pi and append apply assert break class clock complex conjugate continue def del dir dot eigenvectors elif else except exec finally for from global identity if import in is lambda matrixmultiply max mkdir min maximum minimum not or pass path print raise random reduce return remove rename rmdir shape sleep split string system sum time try wait while Data Float Int Numeric array open write close float int input range type zeros " ".
Method,. append method. close method. #!/usr/bin/env python class player: def kick(self): return 'kick the ball' def heading(self): return 'heading the ball' ronaldo=player() print ronaldo.kick() print ronaldo.heading() instantiation
#!/usr/bin/env python class player: strength= 50 speed= 90 stamina= 70 agility= 90 def kick(self): return 'kick the ball' def heading(self): return 'heading the ball' def take_a_rest(self): self.stamina= self.stamina+10 self.agility= self.agility+2 self.speed= self.speed+5 ronaldo=player() print ronaldo.kick() print ronaldo.heading() print 'Ronaldo' print ronaldo.strength print ronaldo.speed print ronaldo.stamina print ronaldo.agility print 'Ronaldo takes a rest' ronaldo.take_a_rest() print 'Ronaldo' print ronaldo.strength print ronaldo.speed print ronaldo.stamina print ronaldo.agility print 'Park' park=player() print park.kick() print park.heading() print park.strength print park.speed print park.stamina print park.agility print 'Park takes a rest' park.take_a_rest() print 'Park' print park.strength print park.speed print park.stamina print park.agility ihlee@cetus0:~/lecture_programming$ class1.py kick the ball heading the ball Ronaldo 50 90 70 90 Ronaldo takes a rest Ronaldo 50 95 80 92 Park kick the ball heading the ball 50 90 70 90 Park takes a rest Park 50 95 80 92 ihlee@cetus0:~/lecture_programming$. method.,.
method #!/usr/bin/env python class Stat: def set_data(self,player,goals,assts): self.player= player self.goals= goals self.assts= assts def printdata(self): print 'player: ', self.player print 'goals: ', self.goals print 'assts: ', self.assts def init (self,player,goals,assts): self.set_data(player,goals,assts) print 'A new Stat class' ron=stat('ronaldo',16,2) ron.set_data('ronaldo',17,3) ron.printdata() lee=stat('lee',0,0) lee.printdata() ihlee@cetus0:~/lecture_programming$ class3.py A new Stat class player: Ronaldo goals: 17 assts: 3 A new Stat class player: Lee goals: 0 assts: 0 ihlee@cetus0:~/lecture_programming$.. self: method. Method,.
#!/usr/bin/env python class Shape: area=0. varia=9. def add (self,other): return self.area + other.area ihlee@cetus0:~/lecture_programming$ class4.py 40.0 40.0 ihlee@cetus0:~/lecture_programming$ a=shape() a.area=10. b=shape() b.area=30. print a+b print a. add (b) method (overload).
#!/usr/bin/env python class Shape: area=0. varia=9. def add (self,other): return self.area + other.area def cmp (self,other): if self.area < other.area : return -1 elif self.area == other.area : return 0 else : return 1 ihlee@cetus0:~/lecture_programming$ class5.py 40.0 40.0-1 b > a ihlee@cetus0:~/lecture_programming$ a=shape() a.area=10. b=shape() b.area=30. print a+b print a. add (b) print a. cmp (b) if a > b : print 'a > b' if a < b : print 'b > a' if a == b : print 'b = a'
#!/usr/bin/env python f=open('abc.txt') buffer=f.read() print buffer f.close() ihlee@cetus0:~/lecture_programming$ file.py I am a boy. You are a girl. ihlee@cetus0:~/lecture_programming$ #!/usr/bin/env python f=open('abc.txt') buffer=f.read() print buffer f.close() g=open('new_abc.txt','w') g.write('my Love') g.close() #!/usr/bin/env python f=open('abc.txt') line=f.readline() print line line=f.readline() print line f.close() ihlee@cetus0:~/lecture_programming$ file3.py I am a boy. You are a girl. ihlee@cetus0:~/lecture_programming$
phi_trj=[[0.]*(iidd) for j in xrange(np+1)] psi_trj=[[0.]*(iidd) for j in xrange(np+1)] for ifile in range(0,np+1): tmp_file=asc2num('rama',ifile)+'.txt' fp_rama=open(tmp_file,'r') jth=0 while 1: line=fp_rama.readline() if not line : break else: token=line.split() phi_trj[ifile][jth]=float(token[0]) psi_trj[ifile][jth]=float(token[1]) jth=jth+1 fp_rama.close() #!/usr/bin/env python from string import atof version ="1.0.0" np=0 f=open('coord.xyz','r') while 1: line=f.readline() if not line : break token=line.split() natom=int(token[0]) line=f.readline() np=np+1 cx= cy= cz=0. for j in range(natom): line=f.readline() token=line.split() cx=cx+float(token[1]) cy=cy+float(token[2]) cz=cz+float(token[3]) cx=cx/float(natom) cy=cy/float(natom) cz=cz/float(natom) format='%21.11f %21.11f %21.11f' print format % (cx,cy,cz) f.close() np=np-1 print natom,np,'natom,np','(',np+1,'frames',')'
#!/usr/bin/env python import os,sys,time,random from posix import system version ="1.0.0" natom=597 nk=0 np=20 amp=5.0e00 pseed=1421 # #natom=int(raw_input("natom :")) #np=int(raw_input("np :")) #nk=int(raw_input("nk :")) #amp=float(raw_input("amp :")) #pseed=int(raw_input("pseed :")) if np <= 0: np=2 if nk <= 0: nk=np-1 print pseed, 'pseed in python' random.seed(pseed) system('grep CA confi.a > calpha_set.dat') ca_site=[] f_ca=open('calpha_set.dat','r') while 1: line=f_ca.readline() if not line : break else: token=line.split() ca_site.append(int(token[0])) f_ca.close() num_calpha=len(ca_site) print num_calpha,' the number of Calpha atoms' print ca_site # akx=[ [[0.]*nk for j in range(3)] for i in range(natom) ] aky=[ [[0.]*nk for j in range(3)] for i in range(natom) ] akz=[ [[0.]*nk for j in range(3)] for i in range(natom) ] for i in range(1,natom+1): for k in range(1,nk+1): akx_t= (random.random()-0.5)*amp/float(k)/float(k) aky_t= (random.random()-0.5)*amp/float(k)/float(k) akz_t= (random.random()-0.5)*amp/float(k)/float(k) akx[i-1][0][k-1]=akx_t aky[i-1][1][k-1]=aky_t akz[i-1][2][k-1]=akz_t # g_sine=open('r_fort.40','w') g_sine.write('%8d %8d %s %s \n' %( nk,natom, 'nk','natom')) istart=0 for i in range(len(ca_site)): ica=ca_site[i] ica=ica-1 if i >0 : istart=ca_site[i-1] print istart+1,ica+1 for iatom in range(istart,ica): for k in range(nk): akx[iatom][0][k]=akx[ica][0][k]+((random.random()-0.5)*1.e-2)/float(k+1)/float(k+1) aky[iatom][1][k]=aky[ica][1][k]+((random.random()-0.5)*1.e-2)/float(k+1)/float(k+1) akz[iatom][2][k]=akz[ica][2][k]+((random.random()-0.5)*1.e-2)/float(k+1)/float(k+1) # ica=ca_site[len(ca_site)-1] print ica+1,natom for iatom in range(ica,natom): for k in range(nk): akx[iatom][0][k]=akx[ica][0][k]+((random.random()-0.5)*1.e-2)/float(k+1)/float(k+1) aky[iatom][1][k]=aky[ica][1][k]+((random.random()-0.5)*1.e-2)/float(k+1)/float(k+1) akz[iatom][2][k]=akz[ica][2][k]+((random.random()-0.5)*1.e-2)/float(k+1)/float(k+1) # for i in range(natom): g_sine.write('\n') for k in range(nk): g_sine.write('%16.8f %16.8f %16.8f \n' %( akx[i][0][k], aky[i][1][k], akz[i][2][k])) g_sine.close()
#!/usr/bin/env python from Numeric import * from string import atof,atoi,split from posix import getpid,rename,remove,mkdir,rmdir,system from shutil import copyfile import time,os,sys,random def prtjobid(): f=open('job_id_printed','w') f.write('%15d \n' % os.getpid()) st= time.asctime() f.write('%23s \n' % st) f.close() return ihlee@cetus0:~/lecture_programming$ time_test.py Thu Jan 17 10:39:36 2008 2003 iseed in python 0.112763723084 0.964931797669 0.693309755783 0.283805968663 0.667882922285 0.790271630375 0.803438096482 0.849369062771 0.707105471223 0.0992155723476 ihlee@cetus0:~/lecture_programming$ print time.asctime() prtjobid() iseed=2003 print iseed, ' iseed in python' random.seed(iseed) for i in range(10): print random.random() time.py. time..
(inheritance) #!/usr/bin/env python class person: eyes=2 nose=1 mouth=1 ears=2 arms=2 legs=2 def eat(self): print 'eat' def sleep(self): print 'sleep' def talk(self): print 'talk' ihlee@cetus0:~/lecture_programming$ class2.py 1 talk 1 talk play football ihlee@cetus0:~/lecture_programming$ class player(person): def football(self): print 'play football' lee=person() print lee.mouth lee.talk() seol=player() print seol.mouth seol.talk() seol.football() Player person. football
,... ` `;? CALL SYSTEM('application.py'). call system. ( ).,.. : / == call system('python_driver.py') ==,. job..
Global variables ihlee@cetus0:~/lecture_programming$ cat glob1.py #!/usr/bin/env python version ="1.0.0" def demo (f_in): global somevar # shared with main code demo.tom = 16 # An attribute accessible from main code somevar += 1 another = 12 # A local variable, independent of main code res = f_in+14 # Value passed in (f_in) return res somevar = 27 # accessed in function via global another = 17 # not accessed in function pval = 16 # accessed in function via parameter print demo(pval),' demo(pval)' print demo.tom,' demo.tom' # function attribute print somevar,'somevar' print another,'another' ihlee@cetus0:~/lecture_programming$ glob1.py 30 demo(pval) 16 demo.tom 28 somevar 17 another ihlee@cetus0:~/lecture_programming$ http://www.wellho.net/forum/programming-in-python-and-ruby/variable-scope-in-python.html
Numeric python; Numpy a powerful N-dimensional array object sophisticated (broadcasting) functions basic linear algebra functions basic Fourier transforms sophisticated random number capabilities tools for integrating Fortran code. tools for integrating C/C++ code. http://www.scipy.org/numpy_example_list_with_doc/, MATLAB numpy.fft numpy.linalg numpy.random numpy.oldnumeric from Numeric import * #When you create an array you must then explicitly tell Python you are doing so as follows: >>> from Numeric import * >>> xx = array([1, 5, 6.5, -11]) >>> print xx [ 1. 5. 6.5-11. ] >>> print xx[0] 1.0 >>> print xx[3] -11.0 >>> xx = zeros(5, Int) >>> print xx [0 0 0 0 0] >>> yy = zeros(4, Float) >>> print yy [ 0. 0. 0. 0.] Numeric -- -- numpy.oldnumeric ihlee@cetus0:~/lecture_programming$ python Python 2.3.5 (#2, Oct 16 2006, 19:59:54) [GCC 3.3.5 (Debian 1:3.3.5-13)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> from Numeric import * >>> m=zeros((3,3)) >>> print m [[0 0 0] [0 0 0] [0 0 0]] >>> print identity(3) [[1 0 0] [0 1 0] [0 0 1]] >>> http://numpy.scipy.org/numpydoc/numpy.html
#!/usr/bin/env python # import os,sys,time,random from Numeric import * from LinearAlgebra import solve_linear_equations from string import atof,atoi,split from posix import getpid,rename,remove,mkdir,rmdir,system from shutil import copyfile version ="1.0.0" def gen_ak_path(np,xvec): # written by In-Ho Lee, KRISS, Jan. 29, 2003. us=0.0 ; uf=1.0 ; du=(uf-us)/float(np) #. bvec=zeros(np-1,float) for j in range(1,np): uu=us+du*float(j) bvec[j-1]=xvec[j]-xvec[0]-uu*(xvec[np]-xvec[0]) amat=zeros((np-1,np-1),float) for j in range(1,np): uu=us+du*float(j) for k in range(1,np): amat[j-1,k-1]=sin(float(k)*pi*uu) # ak_vec=solve_linear_equations(amat,bvec) #! return len(ak_vec),ak_vec
Universal functions >>> from Numeric import * >>> m=zeros((3,3)) >>> print m [[0 0 0] [0 0 0] [0 0 0]] >>> print identity(3) [[1 0 0] [0 1 0] [0 0 1]] >>> sin(m) array([[ 0., 0., 0.], [ 0., 0., 0.], [ 0., 0., 0.]]) >>> m=m+1. >>> print m [[ 1. 1. 1.] [ 1. 1. 1.] [ 1. 1. 1.]] >>> sin(m) array([[ 0.84147098, 0.84147098, 0.84147098], [ 0.84147098, 0.84147098, 0.84147098], [ 0.84147098, 0.84147098, 0.84147098]]) >>> m=m+m >>> print m [[ 2. 2. 2.] [ 2. 2. 2.] [ 2. 2. 2.]] >>> n=m-0.5 >>> print n [[ 1.5 1.5 1.5] [ 1.5 1.5 1.5] [ 1.5 1.5 1.5]] >>> qq=m*n >>> print qq [[ 3. 3. 3.] [ 3. 3. 3.] [ 3. 3. 3.]] >>> pp=matrixmultiply(m,n) >>> print pp [[ 9. 9. 9.] [ 9. 9. 9.] [ 9. 9. 9.]] >>> http://numpy.sourceforge.net/numdoc/html/numdoc.htm#pgfid-57315
#!/usr/bin/env python import os,sys,time from shutil import copyfile from posix import getpid,rename,remove,mkdir,rmdir,system from time import * from Numeric import * from LinearAlgebra import * def test_sort(): n=5 arrin=zeros((n),float) arrin=array([1., 3., 4., 5., 2.]) indx=argsort(arrin) print indx for i in range(5): j=indx[i] print arrin[j] return test_sort() [0 4 1 2 3] 1.0 2.0 3.0 4.0 5.0 #!/usr/bin/env python from Numeric import * from time import * from string import atof,atoi,split from posix import getpid,rename,remove,mkdir,rmdir,system from shutil import copyfile import os,sys,random version ="1.0.0" if name == " main ": nnn=5 at=zeros(nnn,float) s_at=zeros(nnn,float) at[0]=-3. at[1]=3. at[2]=-4. at[3]=-5. at[4]=-6. print at,' before' index=argsort(at) for i in range(nnn): s_at[i]=at[index[i]] print s_at,' after' print index,' index' [-3. 3. -4. -5. -6.] before [-6. -5. -4. -3. 3.] after [4 3 2 0 1] index
Numeric Python (Numpy).,. array1=zeros(100,float) 0.0e0. Numeric-20.3. 100 ( 0,1,2,...99) 0.0e0. (python 0.0 0.0e0. FORTRAN real*8. 0.d0.) double precision. 2 0.0 array2d=zeros((100,100),float). pi. Numeric. " from Numeric import *". "Numeric (*) ". "from LinearAlgebra import solve_linear_equations" solution_vector=solve_linear_equations(a_matrix,b_vector). python. ( ) package. ( C FORTRAN). Numeric python ATLAS(Automatically Tuned Linear Algebra Software; fast linear algebra routines).. ATLAS MAPLE(v7 and higher), MATLAB(v6.0 and higher), Mathematica(forthcoming), MAPLE(v7 and higher), Octave.
vs www.netlib.org C FORTRAN. python. C FORTRAN. python.. (C FORTRAN very-high-level.).., ( ). : : tuple:. ( ) : Numeric python. http://numpy.scipy.org/numpy.pdf
#!/usr/bin/env python #from Numeric import * #from string import atof,atoi,split from string import * #from posix import getpid,rename,remove,mkdir,rmdir,system #from shutil import copyfile #import os,sys,time,random version ="1.0.0" ihlee@cetus0:~/lecture_programming$ gc.py 16.6666666667 ihlee@cetus0:~/lecture_programming$ def gc(dna): return (count(dna, 'c')+count(dna, 'g'))/float(len(dna))*100.0 print gc('atgtaatgatat')
ihlee@cetus0:~/lecture_programming$ python Python 2.3.5 (#2, Oct 16 2006, 19:59:54) [GCC 3.3.5 (Debian 1:3.3.5-13)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> 1+1 2 >>> xx=1+10 >>> print xx 11 >>> yy=xx+xx**2 >>> print yy 132 >>> >>> LL=[1, 2, 3, 4, 5] >>> print LL [1, 2, 3, 4, 5] >>> sum(ll) 15 >>> leng(ll) Traceback (most recent call last): File "<stdin>", line 1, in? NameError: name 'leng' is not defined >>> len(ll) 5 >>> sum(ll)/len(ll) 3 >>> :list append() extend() insert() remove() pop() index() count() sort() reverse()
if os.path.exists('chg') : os.path.exists CHG., os.path.exists('chg') 0 ( ) 'CHG'. system(), system('qsub work.0003/pbs_script_file') pbs. j=3 ; f_name=ascnum('work.',j)+'/pbs_script_file' ; system('qsub '+f_name+'') glob glob.. import glob,time import os.path file_list=glob.glob('*') for f_name in file_list: if os.path.isfile(f_name): print f_name,' is a regular file' print os.path.getsize(f_name) print os.path.getatime(f_name) print time.ctime(os.path.getatime(f_name)) elif os.path.isdir(f_name): print f_name,' is a directory' elif os.path.islink(f_name): print f_name,' is a symbolic link' glob.glob('*.gif') glob.glob('*.dat'). os.listdir('.'). :os.getcwd() www.google.com.. " + python" ("python -snake -monty" ). www.python.org/doc/current/lib. os.getcwd() os.getuid() os.getuname() os.path.getsize(file_name) os.path.split(file_name)
g=open('sad.i','w') # sad.i g.write('%16d %16.8f %s \n' % ( natom,delrr, 'natom,delrr' )) g.write('%16d %16d %16d %s \n' % ( iseed1,iseed2,ktmax, 'iseed1,iseed2,ktmax')) g.write('%16.8f %16.8f %s \n' % ( xs_local,xf_local, 'xs_local,xf_local')) g.write('%16.8f %16.8f %s \n' % ( ys_local,yf_local, 'ys_local,yf_local')) g.write('%16.8f %16.8f %s \n' % ( zs_local,zf_local, 'zs_local,zf_local')) g.close() ihlee@cetus0:~/lecture_programming$ python Python 2.3.5 (#2, Oct 16 2006, 19:59:54) [GCC 3.3.5 (Debian 1:3.3.5-13)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> f=open('input_xx','r') >>> for line in f:... for i in range(len(line.split())):... print line.split()[i]... 1 2. abc >>> ihlee@cetus0:~/lecture_programming$ vi input_xx ihlee@cetus0:~/lecture_programming$ cat input_xx 1 2. abc ihlee@cetus0:~/lecture_programming$.close().flush().write() d e f f (e, f) c format='%16.8e %16.8e %16.8e %16.8e %12.5e %12.5e %22s' print format %(ac+cc+zz,ac,cc,zz,onma,etar,'o,a,p1,p2,onma,etarget')
atof, atoi ihlee@cetus0:~/lecture_programming$ cat input_xx 1 2. abc ihlee@cetus0:~/lecture_programming$ cat read_test.py #!/usr/bin/env python from string import atof,atoi,split f=open("input_xx",'r') for line in f: int_aaa=atoi(line.split()[0]) rel_bbb=atof(line.split()[1]) chr_ccc= line.split()[2] print int_aaa print rel_bbb print chr_ccc f.close() ihlee@cetus0:~/lecture_programming$ read_test.py 1 2.0 abc ihlee@cetus0:~/lecture_programming$ int, float, complex raw_input().
",, midpoint" def read_energy(fname): # written by In-Ho Lee, KRISS, Jan. 12, 2002. f=open(fname,'r') for line in f: for i in range(len(line.split())): if line.split()[i] == 'midpoint': energy= atof(line.split()[i-2]) f.close() return energy python 2.0. def read_energy_saddle(fname): # written by In-Ho Lee, KRISS, March 6, 2003. select=1 f=open(fname,'r') while select : line=f.readline() for i in range(len(line.split())): if line.split()[i] == 'midpoint': energy= atof(line.split()[i-2]) select=0 break f.close() return energy
PBS (: PBS. wwait. PBS CPU.) system('qsub ursa_vasp') time.sleep(10.) wwait('done_vasp') PBS job touch DONE_VASP DONE_VASP python ( ). (touch. 0.) python wwait DONE_VASP ( 60.).,. job, job serial job. system('qsub work.0003/pbs_script_file') pbs. j=3 ; f_name=ascnum('work.',j)+'/pbs_script_file' ; system('qsub '+f_name+'')
#!/usr/bin/env python from Numeric import * # from string import atof,atoi,split from posix import getpid,rename,remove,mkdir,rmdir,system from shutil import copyfile import os,sys,time,random version ="1.0.0" def prtjobid(): # written by In-Ho Lee, KRISS, Jan. 12, 2002. f=open('job_id_printed','w') f.write('%15d \n' % os.getpid()) # st=time.asctime() # f.write('%23s \n' % st) f.close() return 6291 Tue Jan 15 11:48:19 2008,. ps u ihlee, process ID kill 9 process ID. prtjobid() ihlee@cetus0:~/lecture_programming$ ls -ltra.
Standard input #!/usr/bin/env python #from Numeric import * #from string import atof,atoi,split #from posix import getpid,rename,remove,mkdir,rmdir,system #from shutil import copyfile #import os,sys,time,random version ="1.0.0" print "Please give a number: " a = input() print "And another: " b = input() print "The sum of these numbers is: " print a + b ihlee@cetus0:~/lecture_programming$ input.py Please give a number: 1 And another: 2 The sum of these numbers is: 3 ihlee@cetus0:~/lecture_programming$ raw_input() sys.exit() #!/usr/bin/env python import sys print "Type a value: ", # comma prevents newline value = sys.stdin.readline() # use stdin explicitly print value sys.stdout.write("hello world\n") # \n= newline ihlee@cetus0:~/lecture_programming$ stdin.py Type a value: 3 3 ihlee@cetus0:~/lecture_programming$ stdin.py Type a value: 33. 33. ihlee@cetus0:~/lecture_programming$ stdin.py Type a value: 'abc' 'abc' ihlee@cetus0:~/lecture_programming$ stdin.py Type a value: "defg" "defg"
>>> import sys >>> dir(sys) [' displayhook ', ' doc ', ' excepthook ', ' name ', ' stderr ', ' stdin ', ' stdout ', '_getframe', 'api_version', 'argv', 'builtin_modul e_names', 'byteorder', 'call_tracing', 'callstats', 'copyright', 'displayhook', 'exc_clear', 'exc_info', 'exc_type', 'excepthook', 'exec_prefix', 'executable', 'exit', 'getcheckinterval', 'getdefaultencoding', 'getdlopenflags', 'getfilesystemencoding', 'getrecursionlimit', 'getrefcount', 'hexversion', 'last_tracebac k', 'last_type', 'last_value', 'maxint', 'maxunicode', 'meta_path', 'modules', 'path', 'path_hooks', 'path_importer_cache', 'platform', 'prefix', 'ps1', 'ps 2', 'setcheckinterval', 'setdlopenflags', 'setprofile', 'setrecursionlimit', 'settrace', 'stderr', 'stdin', 'stdout', 'version', 'version_info', 'warnoptions'] >>>
>>> x = 5 >>> print x 5 >>> a = b = c = 0 >>> myfloat = 0.1 >>> myfloat = 2.0 >>> myfloat = 3.14159256 >>> myfloat = 1.6e-19 >>> myfloat = 3e8 >>> print range(0, 6) [0, 1, 2, 3, 4, 5] >>> print range(10, 20, 2) [10, 12, 14, 16, 18] >>> print a, b, c 0 0 0 >>> mystring = "Here is a string" >>> mystring = 'Here is another' >>> i = 2 >>> i = i + 1 >>> print i 3 >>> print (2 * (3-1)) * 4 16 >>> myinteger = 0 >>> myinteger = 15 >>> myinteger = -23 >>> myinteger = 2378 sumsquares = 0 # sumsquares must have a value because we increment # it later. for i in [0, 1, 2, 3, 4, 5]: print "i now equal to:", i sumsquares = sumsquares + i**2 # sumsquares incremented here print "sum of squares now equal to:", sumsquares print "------" print "Done."
,. FORTRAN 90 vector matrix. >>> range(1,10,2) [1, 3, 5, 7, 9] >>> range(10,1,-1) [10, 9, 8, 7, 6, 5, 4, 3, 2] >>> range(10,0,-1) [10, 9, 8, 7, 6, 5, 4, 3, 2, 1] >>> >>> istart=1; ifinish=9 ; idirection=1 >>> for i in range(istart,ifinish+1,idirection):... print i... 1 2 3 4 5 6 7 8 9 >>> istart=9; ifinish=-1; idirection=-1 >>> for i in range(istart,ifinish+1,idirection):... print i... 9 8 7 6 5 4 3 2 1 >>>
from Numeric import * >>> a=array((1,2)) >>> b=array((3,4)) >>> dot(a,b) 11 >>> a=array(((1,2),(3,4))) >>> b=a >>> matrixmultiply(a,b) #! array([[ 7, 10], [15, 22]]) >>> sin(complex(1,3.)) (8.4716454543001483+5.4126809231781934j) >>> >>> for i in [1,2,3,4]:... print i... 1 2 3 4.. http://numpy.scipy.org/numpy.pdf
>>> z=complex(1.,2.) >>> z.conjugate() (1-2j) >>> z.real 1.0 >>> z.imag 2.0 >>> arange(0.,1.,0.20) array([ 0., 0.2, 0.4, 0.6, 0.8]) >>> arange(0.0,1.2,0.20) array([ 0., 0.2, 0.4, 0.6, 0.8, 1. ]) >>> >>> q=arange(0.0,1.0,0.2) >>> sin(q) array([ 0., 0.19866933, 0.38941834, 0.56464247, 0.71735609]) >>> >>> max(0.,2.,3.) 3.0 >>> maximum(0.,2.) 2.0 >>> min(0.,-1,-3.) -3.0 >>> minimum(-3.,3.) -3.0 >>> a=[1,2,3.,4] >>> sum(a) 10.0
>>> b=[4,3,7,9] >>> sort(b) array([3, 4, 7, 9]) >>> amat=zeros((2,2),float) >>> amat[0][0]=1 >>> amat[1][1]=1 >>> trace(amat) 2.0 >>> from LinearAlgebra import * >>> determinant(amat) 1.0 >>> eigenvalues(amat) array([ 1., 1.]) >>> eigenvectors(amat) (array([ 1., 1.]), array([[ 1., 0.], [ 0., 1.]])) >>> b=inverse(amat) >>> print b [[ 1. 0.] [ 0. 1.]] >>> >>> b=zeros(2,float) >>> solve_linear_equations(amat,b) array([ 0., 0.]) >>> >>> transpose(amat) array([[ 1., 0.], [ 0., 1.]]) >>> >>> argmax(b) 0 >>> argmin(b) 0 #.. List :list Seq.append() Seq.extend() Seq.insert() Seq.remove() Seq.pop() Seq.index() Seq.count() Seq.sort() Seq.reverse() append method. Sequences Strings (type str) of length 0, 1, 2 '', '1', "12", 'hello\n' Tuples (type tuple) of length 0, 1, 2, etc: () (1,) (1,2) # parentheses are optional if len > 0 Lists (type list) of length 0, 1, 2, etc: [] [1] [1,2]
#!/usr/bin/env python #from Numeric import * #from string import atof,atoi,split #from posix import getpid,rename,remove,mkdir,rmdir,system #from shutil import copyfile #import os,sys,time,random version ="1.0.0" sumsquares = 0 # sumsquares must have a value because we increment # it later. for i in [0, 1, 2, 3, 4, 5]: print "i now equal to:", i sumsquares = sumsquares + i**2 # sumsquares incremented here print "sum of squares now equal to:", sumsquares print "------" print "Done." i now equal to: 0 sum of squares now equal to: 0 ------ i now equal to: 1 sum of squares now equal to: 1 ------ i now equal to: 2 sum of squares now equal to: 5 ------ i now equal to: 3 sum of squares now equal to: 14 ------ i now equal to: 4 sum of squares now equal to: 30 ------ i now equal to: 5 sum of squares now equal to: 55 ------ Done. ihlee@cetus0:~/lecture_programming$
(matrix) (diagonalization). H 2 (n by n). H=T+V # Hamiltonian setup val,vec=eigenvectors(h) # matrix diagonalization (eigenvalue/vector!) val,vec=ev_sort(val,vec) # solution sorting def ev_sort(eigval,eigvec): newval=zeros(len(eigval),float) newvec= zeros(eigvec.shape,float) index=argsort(eigval) for i in index: newval[i]=eigval[index[i]] newvec[i,:]=eigvec[index[i],:] # return newval,newvec. C, Fortran.,.. aaa[i,j] aaa[i][j]. Numpy [i,j], list, [i][j].
>>> def addnumbers(x, y): sum = x + y return sum >>> x = addnumbers(5, 10) >>> print x 15 >>> def addnumbers(x, y): sum = x + y x = 100000 return [sum, x] >>> x = 5 >>> y = 10 >>> answer = addnumbers(x, y) >>> print answer[0] 15 >>> print answer[1] 100000 iseed=2003 print iseed, ' iseed in python' random.seed(iseed) for i in range(10): print random.random() # # 10 for i in range(2000): # 2000 gen_sad_input() # sad.i system('sad.x >TMPFILE') # sad.x copyfile('dim.dat', ascnum('dimer.',i)) copyfile('mid.xyz', ascnum_xyz('midpt.',i)) # system('rm dim.dat') # system('rm TMPFILE')
>>> fin = open("input.dat", "r") >>> line = fin.readline() >>> print line 1 5.06 78 15 def gen_sad_input(): natom=60 # ktmax=2000 delrr=1.e-1 xs_local= 4. -random.random()*2. xf_local= 8. +random.random()*2. ys_local=-1. -random.random()*2. yf_local= 1. +random.random()*2. zs_local=-3. -random.random()*2. zf_local= 3. +random.random()*2. iseed1= int( (random.random())*31327.0e0+0. ) iseed2 = int( (random.random())*30080.0e0+0. ) g=open('sad.i','w') # sad.i g.write('%16d %16.8f %s \n' % ( natom,delrr, 'natom,delrr' )) g.write('%16d %16d %16d %s \n' % ( iseed1,iseed2,ktmax, 'iseed1,iseed2,ktmax')) g.write('%16.8f %16.8f %s \n' % ( xs_local,xf_local, 'xs_local,xf_local')) g.write('%16.8f %16.8f %s \n' % ( ys_local,yf_local, 'ys_local,yf_local')) g.write('%16.8f %16.8f %s \n' % ( zs_local,zf_local, 'zs_local,zf_local')) g.close() return >>> data = split(line) >>> print data ['1', '5.06', '78', '15'] >>> x = int(data[0]) >>> print x 1 >>> y = float(data[1]) >>> print y 5.06
ihlee@cetus0:~/lecture_programming$ python Python 2.3.5 (#2, Oct 16 2006, 19:59:54) [GCC 3.3.5 (Debian 1:3.3.5-13)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> abc='abc' >>> print abc[0] a >>> print abc[1] b >>> print abc[2] c >>> princ abc[3] File "<stdin>", line 1 princ abc[3] ^ SyntaxError: invalid syntax >>> print abc[-1] c >>> print abc[-2] b >>> print abc[-3] a >>> print abc[-4] Traceback (most recent call last): File "<stdin>", line 1, in? IndexError: string index out of range >>> len(abc) 3 >>> print sentence1[2:] is is a book. >>> >>> sentence1='this is a book.' >>> print sentence1 This is a book. >>> print sentence1[2:10] is is a >>> print sentence1[2:12] is is a bo >>> len(sentence1) 15,,, tuple
#!/usr/bin/env python from string import * def how_many_words(s): list = split(s) return len(list) s='energy = -100.000 ev' print how_many_words(s) 4 def find(str, ch): index = 0 while index < len(str): if str[index] == ch: return index index = index + 1 return -1 >>> astring='abcdefg' >>> bstring=astring[:-1] >>> print bstring abcdef >>> print astring[:] abcdefg >>> print astring[1:] bcdefg >>> print astring[3:] defg
Join and Split #!/usr/bin/env python days = ["Monday","Tuesday","Webnesday"] all = "<br>".join(days) print all print ''.join(days) if "o" in "Melksham": print "yes" ihlee@cetus0:~/lecture_programming$ join_split.py Monday<br>Tuesday<br>Webnesday MondayTuesdayWebnesday ['This', 'is', 'a', 'line', 'of', 'text'] Smiota Tom ihlee@cetus0:~/lecture_programming$ say = "This is a line of text" parts = say.split(" ") first, second = "Tom Smiota".split(" ") second, first = first, second print parts print first print second ( ), join.
s.capitalize() s.count() s.endswith() s.strip() s.split() s.swapcase() s.upper() Python 2.4 Quick Reference Python reference card http://www.digilife.be/quickreferences/qrc/python%202.4%20quick%20reference%20card.pdf http://rgruet.free.fr/pqr2.3.html
#!/usr/bin/env python from Numeric import * from string import atof,atoi,split from posix import getpid,rename,remove,mkdir,rmdir,system from shutil import copyfile import os,sys,time,random version ="1.0.0" if name == " main ": # # *.bin # convert.x #. convert.x "convert.x test.bin test". # # yeskim98, March/26/2003 # system('ls *.bin >namelist') # f=open('namelist','r') for line in f: name=line.split()[0] name1=name[:-4] # 4 system('./convert.x '+name+' '+name1+'') f.close() abc[:] : abc[::-1]
n1 n 1 f(n). f(13)=6. f(n)=n 1.? #!/usr/bin/env python #from Numeric import * #from string import atof,atoi,split #from posix import getpid,rename,remove,mkdir,rmdir,system #from shutil import copyfile #import os,sys,time,random version ="1.0.0" counter=0 acc = 0 for n in xrange(1,10000000000): num1 = str(n).count('1') if num1 > 0 : acc = acc+num1 if acc == n: counter=counter+1 print "counter = ",counter,n if counter == 2 : break ihlee@cetus0:~/lecture_programming$ google0.py counter = 1 1 counter = 2 199981 ihlee@cetus0:~/lecture_programming$ one by one, passing each to the variable n in turn range vs xrange >>> range(1,10) [1, 2, 3, 4, 5, 6, 7, 8, 9] >>> xrange(1,10) xrange(1, 10) >>> http://incredible.egloos.com/2980264
num1 = str(n).count('1') num1 = str(n) length=len(num1) num1num=0 for jj in range(length): if num1[jj] == '1': num1num=num1num+1 num1=num1num
{dictionary} #!/usr/bin/env python author = {"php":"rasmus Lerdorf",\ "perl":"larry Wall",\ "python":"guido van Rossum"} ihlee@cetus0:~/lecture_programming$ dic1.py Python was written by Guido van Rossum php was written by Rasmus Lerdorf perl was written by Larry Wall ihlee@cetus0:~/lecture_programming$ print "Python was written by",author["python"] print "php was written by",author["php"] print "perl was written by",author["perl"] #!/usr/bin/env python author = {"php":"rasmus Lerdorf",\ "perl":"larry Wall",\ "tcl":"john Ousterhout",\ "awk":"brian Kernighan",\ "java":"sun Microsystems",\ "parrot":"simon Cozens",\ "python":"guido van Rossum"} author["c"] = "Brian Kernighan" del author["parrot"] author["java"] = "James Gosling" ihlee@cetus0:~/lecture_programming$ dic2.py Python was written by Guido van Rossum c is the child of Brian Kernighan java is the child of James Gosling python is the child of Guido van Rossum perl is the child of Larry Wall tcl is the child of John Ousterhout awk is the child of Brian Kernighan php is the child of Rasmus Lerdorf ihlee@cetus0:~/lecture_programming$ print "Python was written by",author["python"] print for language in author.keys(): print language,"is the child of",author[language] http://www.wellho.net/course/python.html
Tuple,.. immutable Tuple. >>> tpl1=(1,2) >>> tpl2=('a', 'b') >>> print tpl1 (1, 2) >>> print tpl2 ('a', 'b') >>> tpl3=tpl1+tpl2 >>> print tpl3 (1, 2, 'a', 'b') >>> >>> for i in abc:... print i... a b c >>> for i in positions:... print i... GK DF MF FW >>> for i in tpl3:... print i... 1 2 a b >>> range(0,10) [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] >>> range(10) [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] >>> for i in ss:... print i... GK DF MF FW ronaldo park rooney lee,, :..
Tuple #!/usr/bin/env python demo_tpl = (1, 3, 6, 10, 15, 21, "lots") for i in range(len(demo_tpl)): print i,"--->",demo_tpl[i] ihlee@cetus0:~/lecture_programming$ tpl1.py 0 ---> 1 1 ---> 3 2 ---> 6 3 ---> 10 4 ---> 15 5 ---> 21 6 ---> lots ihlee@cetus0:~/lecture_programming$ #!/usr/bin/env python values = (100,200,300,500,1000,\ 2000,4000,8000,16000,32000,\ 64000,125000,250000,500000,\ 1000000) being_values = [100,200,300,500,1000,\ 2000,4000,8000,16000,32000,\ 64000,125000,250000,500000,\ 1000000] if values == being_values: print "One" done = tuple(being_values) if values == done: print "Two" Two >>> tu=tuple(range(1,4)) >>> print tu (1, 2, 3) >>> ple=tuple(range(4,7)) >>> tuple1=tu+ple >>> print tuple1 (1, 2, 3, 4, 5, 6) >>> print tuple1[1] 2 >>> tuple2=tu,ple >>> print tuple2 ((1, 2, 3), (4, 5, 6)) >>> print tuple2[1][1] 5 >>>
Tuple #!/usr/bin/env python monthno = [] dayno = [] mlen = (31,28,31,30,31,30,31,31,30,31,30,31) for currentmonth in range (len(mlen)): for day in range (mlen[currentmonth]): dayno.append(day+1) monthno.append(currentmonth+1) # dayno += [day+1] # monthno += [currentmonth+1] wantday = int(raw_input("which day number in year? ")) - 1 print "That's day",dayno[wantday],"of",monthno[wantday] >>> ttt=tuple(range(1,6)) >>> b1,b2,b3,b5 = ttt Traceback (most recent call last): File "<stdin>", line 1, in? ValueError: unpack tuple of wrong size >>> b1,b2,b3,b4,b5 = ttt >>> print ttt (1, 2, 3, 4, 5) >>> print b1,b2,b3,b4,b5 1 2 3 4 5 >>> ihlee@cetus0:~/lecture_programming$ tup_month.py Which day number in year? 31 That's day 31 of 1 ihlee@cetus0:~/lecture_programming$ tup_month.py Which day number in year? 365 That's day 31 of 12 ihlee@cetus0:~/lecture_programming$ tup_month.py Which day number in year? 33 That's day 2 of 2 ihlee@cetus0:~/lecture_programming$ >>> aaa, bbb = -1, 1 >>> print aaa,bbb -1 1 >>> Python 2.3.5 (#2, Oct 16 2006, 19:59:54) [GCC 3.3.5 (Debian 1:3.3.5-13)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> a=[] >>> a+=[0] >>> print a [0] >>> a+=[11] >>> print a [0, 11] >>> a.append(22) >>> print a [0, 11, 22] >>> http://www.wellho.net/course/python.html
,,.. >>> names=['ronaldo', 'baggio', 'rooney'] >>> print names ['ronaldo', 'baggio', 'rooney'] >>> print names[0] ronaldo >>> print names[1:] ['baggio', 'rooney'] >>> names[1]='park' >>> print names ['ronaldo', 'park', 'rooney'] >>> names.append('lee') >>> print names ['ronaldo', 'park', 'rooney', 'lee'] >>> len(names) 4 >>> positions=['gk', 'DF', 'MF', 'FW'] >>> ss=positions+names >>> print ss ['GK', 'DF', 'MF', 'FW', 'ronaldo', 'park', 'rooney', 'lee'] >>> >>> tt=[1, 2, 3] >>> tt3=3*tt >>> print tt3 [1, 2, 3, 1, 2, 3, 1, 2, 3] >>> qq=tt3+names >>> print qq [1, 2, 3, 1, 2, 3, 1, 2, 3, 'ronaldo', 'park', 'rooney', 'lee']
An example that uses most of the list methods: ihlee@cetus0:~/lecture_programming$ python Python 2.3.5 (#2, Oct 16 2006, 19:59:54) [GCC 3.3.5 (Debian 1:3.3.5-13)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> alist=[-1, 3, -7, 9, 8, 4, 4] >>> print alist.count(3) 1 >>> print alist.count(4) 2 >>> print alist.count(11) 0 >>> alist.insert(1,-9) >>> alist [-1, -9, 3, -7, 9, 8, 4, 4] >>> alist.append(88) >>> alist [-1, -9, 3, -7, 9, 8, 4, 4, 88] >>> alist.index(8) 5 >>> alist.remove(88) >>> alist [-1, -9, 3, -7, 9, 8, 4, 4] >>> alist.reverse() >>> alist [4, 4, 8, 9, -7, 3, -9, -1] >>> alist.sort() >>> alist [-9, -7, -1, 3, 4, 4, 8, 9] >>> :list append() extend() insert() remove() pop() index() count() sort() reverse()
String, List, Dictionary, and Class : (1/2) #!/usr/bin/env python class DNA: """Class representing DNA as a string sequence.""" basecomplement = {'A': 'T', 'C': 'G', 'T': 'A', 'G': 'C'} def init (self, s): """Create DNA instance initialized to string s.""" self.seq = s def transcribe(self): """Return as rna string.""" return self.seq.replace('t', 'U') String replace List reverse def reverse(self): """Return dna string in reverse order.""" letters = list(self.seq) letters.reverse() return ''.join(letters) def complement(self): """Return the complementary dna string.""" letters = list(self.seq) letters = [self.basecomplement[base] for base in letters] return ''.join(letters) http://www.onlamp.com/pub/a/python/2002/10/17/biopython.html?page=5
String, List, Dictionary, and Class : (2/2) def reversecomplement(self): """Return the reverse complement of the dna string.""" letters = list(self.seq) letters.reverse() letters = [self.basecomplement[base] for base in letters] return ''.join(letters) def gc(self): """Return the percentage of dna composed of G+C.""" s = self.seq gc = s.count('g') + s.count('c') return gc * 100.0 / len(s) String replace List reverse def codons(self): """Return list of codons for the dna string.""" s = self.seq end = len(s) - (len(s) % 3) - 1 codons = [s[i:i+3] for i in range(0, end, 3)] return codons http://www.onlamp.com/pub/a/python/2002/10/17/biopython.html?page=5
#!/usr/bin/env python class DNA: """Class representing DNA as a string sequence.""" basecomplement = {'A': 'T', 'C': 'G', 'T': 'A', 'G': 'C'} def init (self, s): """Create DNA instance initialized to string s.""" self.seq = s String, List, Dictionary, and Class : def transcribe(self): """Return as rna string.""" return self.seq.replace('t', 'U') def reverse(self): """Return dna string in reverse order.""" letters = list(self.seq) letters.reverse() return ''.join(letters) String replace List reverse def complement(self): """Return the complementary dna string.""" letters = list(self.seq) letters = [self.basecomplement[base] for base in letters] return ''.join(letters) def reversecomplement(self): """Return the reverse complement of the dna string.""" letters = list(self.seq) letters.reverse() letters = [self.basecomplement[base] for base in letters] return ''.join(letters) def gc(self): """Return the percentage of dna composed of G+C.""" s = self.seq gc = s.count('g') + s.count('c') return gc * 100.0 / len(s) def codons(self): """Return list of codons for the dna string.""" s = self.seq end = len(s) - (len(s) % 3) - 1 codons = [s[i:i+3] for i in range(0, end, 3)] return codons http://www.onlamp.com/pub/a/python/2002/10/17/biopython.html?page=5
String, List, Dictionary, and Class : ihlee@cetus0:~/lecture_programming$ vi dna_test.py ihlee@cetus0:~/lecture_programming$ python Python 2.3.5 (#2, Oct 16 2006, 19:59:54) [GCC 3.3.5 (Debian 1:3.3.5-13)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> from DNA import * Traceback (most recent call last): File "<stdin>", line 1, in? ImportError: No module named DNA >>> from dna_test import * >>> dna1=dna('cgacaatac') >>> dna1.transcribe() 'CGACAAUAC' >>> dna1.reverse() 'CATAACAGC' >>> dna1.complement() 'GCTGTTATG' >>> dna1.reversecomplement() 'GTATTGTCG' >>> dna1.gc() 44.444444444444443 >>> dna1.codons() ['CGA', 'CAA', 'TAC'] >>>
from import. from Numeric import sin, cos, tan, sqrt int(),round(),floor(),ceil(),,,,. floor ceil math. float(). dir(module_name)., - -. :,,. A.B A B.
ihlee@cetus0:~/lecture_programming$ python Python 2.3.5 (#2, Oct 16 2006, 19:59:54) [GCC 3.3.5 (Debian 1:3.3.5-13)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> import fibo : from fibo import * >>> a=10 >>> b=1000 fibo.py in lecture_programming directory >>> fib(1000) Traceback (most recent call last): File "<stdin>", line 1, in? NameError: name 'fib' is not defined >>> fibo.fib(1000) 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 >>> print a 10 >>> print b 1000 >>> abc=fibo.fib2(1000) >>> print abc [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987] >>> print a 10 >>> print b 1000 >>>
..... (, )..... (.) Python Library Reference http://docs.python.org/lib/lib.html
#!/usr/bin/env python def fib(n): a, b = 0, 1 while b < n: print b, a, b = b, a+b def fib2(n): result=[] a, b = 0, 1 while b < n : result.append(b) # print b, a, b = b, a+b return result ihlee@cetus0:~/lecture_programming$ python Python 2.3.5 (#2, Oct 16 2006, 19:59:54) [GCC 3.3.5 (Debian 1:3.3.5-13)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> import fibo >>> fibo.fib(1000) 1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 >>> fibo.fib2(1000) [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610, 987] >>> #n=100 #fib(n) #print #n=100 #abc=fib2(n) #print abc fibo.py in lecture_programming directory,.
ihlee@cetus0:~/lecture_programming$ python Python 2.3.5 (#2, Oct 16 2006, 19:59:54) [GCC 3.3.5 (Debian 1:3.3.5-13)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> from fibo import * >>> fib(500) 1 1 2 3 5 8 13 21 34 55 89 144 233 377 >>> fibo.py in lecture_programming directory...py.
BLOSUM50 5-2 -1-2 -1-1 -1 0-2 -1-2 -1-1 -3-1 1 0-3 -2 0-2 -1-1 -2 7-1 -2-4 1 0-3 0-4 -3 3-2 -3-3 -1-1 -3-1 -3-1 0-1 -1-1 7 2-2 0 0 0 1-3 -4 0-2 -4-2 1 0-4 -2-3 4 0-1 -2-2 2 8-4 0 2-1 -1-4 -4-1 -4-5 -1 0-1 -5-3 -4 5 1-1 -1-4 -2-4 13-3 -3-3 -3-2 -2-3 -2-2 -4-1 -1-5 -3-1 -3-3 -2-1 1 0 0-3 7 2-2 1-3 -2 2 0-4 -1 0-1 -1-1 -3 0 4-1 -1 0 0 2-3 2 6-3 0-4 -3 1-2 -3-1 -1-1 -3-2 -3 1 5-1 0-3 0-1 -3-2 -3 8-2 -4-4 -2-3 -4-2 0-2 -3-3 -4-1 -2-2 -2 0 1-1 -3 1 0-2 10-4 -3 0-1 -1-2 -1-2 -3 2-4 0 0-1 -1-4 -3-4 -2-3 -4-4 -4 5 2-3 2 0-3 -3-1 -3-1 4-4 -3-1 -2-3 -4-4 -2-2 -3-4 -3 2 5-3 3 1-4 -3-1 -2-1 1-4 -3-1 -1 3 0-1 -3 2 1-2 0-3 -3 6-2 -4-1 0-1 -3-2 -3 0 1-1 -1-2 -2-4 -2 0-2 -3-1 2 3-2 7 0-3 -2-1 -1 0 1-3 -1-1 -3-3 -4-5 -2-4 -3-4 -1 0 1-4 0 8-4 -3-2 1 4-1 -4-4 -2-1 -3-2 -1-4 -1-1 -2-2 -3-4 -1-3 -4 10-1 -1-4 -3-3 -2-1 -2 1-1 1 0-1 0-1 0-1 -3-3 0-2 -3-1 5 2-4 -2-2 0 0-1 0-1 0-1 -1-1 -1-2 -2-1 -1-1 -1-2 -1 2 5-3 -2 0 0-1 0-3 -3-4 -5-5 -1-3 -3-3 -3-2 -3-1 1-4 -4-3 15 2-3 -5-2 -3-2 -1-2 -3-3 -1-2 -3 2-1 -1-2 0 4-3 -2-2 2 8-1 -3-2 -1 0-3 -3-4 -1-3 -3-4 -4 4 1-3 1-1 -3-2 0-3 -1 5-4 -3-1 -2-1 4 5-3 0 1-1 0-4 -4 0-3 -4-2 0 0-5 -3-4 5 2-1 -1 0 0 1-3 4 5-2 0-3 -3 1-1 -4-1 0-1 -2-2 -3 2 5-1 -1-1 -1-1 -2-1 -1-2 -1-1 -1-1 -1-2 -2-1 0-3 -1-1 -1-1 -1
BLOSUM62 4-1 -2-2 0-1 -1 0-2 -1-1 -1-1 -2-1 1 0-3 -2 0-2 -1 0-1 5 0-2 -3 1 0-2 0-3 -2 2-1 -3-2 -1-1 -3-2 -3-1 0-1 -2 0 6 1-3 0 0 0 1-3 -3 0-2 -3-2 1 0-4 -2-3 3 0-1 -2-2 1 6-3 0 2-1 -1-3 -4-1 -3-3 -1 0-1 -4-3 -3 4 1-1 0-3 -3-3 9-3 -4-3 -3-1 -1-3 -1-2 -3-1 -1-2 -2-1 -3-3 -2-1 1 0 0-3 5 2-2 0-3 -2 1 0-3 -1 0-1 -2-1 -2 0 3-1 -1 0 0 2-4 2 5-2 0-3 -3 1-2 -3-1 0-1 -3-2 -2 1 4-1 0-2 0-1 -3-2 -2 6-2 -4-4 -2-3 -3-2 0-2 -2-3 -3-1 -2-1 -2 0 1-1 -3 0 0-2 8-3 -3-1 -2-1 -2-1 -2-2 2-3 0 0-1 -1-3 -3-3 -1-3 -3-4 -3 4 2-3 1 0-3 -2-1 -3-1 3-3 -3-1 -1-2 -3-4 -1-2 -3-4 -3 2 4-2 2 0-3 -2-1 -2-1 1-4 -3-1 -1 2 0-1 -3 1 1-2 -1-3 -2 5-1 -3-1 0-1 -3-2 -2 0 1-1 -1-1 -2-3 -1 0-2 -3-2 1 2-1 5 0-2 -1-1 -1-1 1-3 -1-1 -2-3 -3-3 -2-3 -3-3 -1 0 0-3 0 6-4 -2-2 1 3-1 -3-3 -1-1 -2-2 -1-3 -1-1 -2-2 -3-3 -1-2 -4 7-1 -1-4 -3-2 -2-1 -2 1-1 1 0-1 0 0 0-1 -2-2 0-1 -2-1 4 1-3 -2-2 0 0 0 0-1 0-1 -1-1 -1-2 -2-1 -1-1 -1-2 -1 1 5-2 -2 0-1 -1 0-3 -3-4 -4-2 -2-3 -2-2 -3-2 -3-1 1-4 -3-2 11 2-3 -4-3 -2-2 -2-2 -3-2 -1-2 -3 2-1 -1-2 -1 3-3 -2-2 2 7-1 -3-2 -1 0-3 -3-3 -1-2 -2-3 -3 3 1-2 1-1 -2-2 0-3 -1 4-3 -2-1 -2-1 3 4-3 0 1-1 0-3 -4 0-3 -3-2 0-1 -4-3 -3 4 1-1 -1 0 0 1-3 3 4-2 0-3 -3 1-1 -3-1 0-1 -3-2 -2 1 4-1 0-1 -1-1 -2-1 -1-1 -1-1 -1-1 -1-1 -2 0 0-2 -1-1 -1-1 -1
Dynamic programming (Python) 1/4 #!/usr/bin/python # Needleman - Wunsch dynamic programming # programmed by Keehyoung Joo import sys, string, os import Numeric ############################################################## MTX = 'BLOSUM62' #MTX = 'PAM70' d = 10 # open gap penalty Amino = ['A','R','N','D','C','Q','E','G','H','I','L', 'K','M','F','P','S','T','W','Y','V','B','Z','X'] # J, O, U not in Amino list ############################################################## def read_score_matrix(s): f=open(mtx,'r') for x in Amino: S[x]={} line = f.readline() values = map(int,string.split(line)) for i in range(len(amino)): y = Amino[i] S[x][y] = values[i]
def align(p,q,v): P = [] Q = [] i = len(p) j = len(q) while 1: if i==j==0: break if V[i,j] == 0: P.insert(0,p[i-1]) Q.insert(0,q[j-1]) i, j = i-1, j-1 elif V[i,j] == 1: P.insert(0,p[i-1]) Q.insert(0,'-') i = i-1 else: P.insert(0,'-') Q.insert(0,q[j-1]) j = j-1 P = "".join(p) Q = "".join(q) return P, Q Dynamic programming (Python) 2/4
Dynamic programming (Python) 3/4 def Dynamic(p,q,Score): # Needleman - Wunsch Algorithm n=len(p) m=len(q) F=Numeric.zeros((n+1,m+1),Numeric.Int) V=Numeric.zeros((n+1,m+1),Numeric.Int) for i in range(1,n+1): F[i,0] = -i * d # d: open gap penalty V[i,0] = 1 for j in range(1,m+1): F[0,j] = -j * d V[0,j] = 2 for i in range(1,n+1): for j in range(1,m+1): L = [F[i-1,j-1]+S[p[i-1]][q[j-1]], F[i-1,j] - d, F[i,j-1] - d] F[i,j] = max(l) V[i,j] = L.index(F[i,j]) P, Q = align(p,q,v) return P, Q, F[n,m]
Dynamic programming (Python) 4/4 ############################################################## if name == ' main ': S={} read_score_matrix(s) p=raw_input('input your seq 1 : ') q=raw_input('input your seq 2 : ') p=p.upper() q=q.upper() Score = 0 (P, Q, Score) = Dynamic(p,q,Score) print print 'Alignment result...\n' print 'Score = %d' % Score print 'seq 1 : %s' % P print 'seq 2 : %s' % Q print
Ramachandran map (Python) :p. 175 # Homework #5 Ramachandran_Map # 2007. 04.11 by Shin, Seung-Woo # Import PROSS module from PROSS import Protein from PROSS import read_pdb # Read PDB file #f_name_pdb=raw_input("input PDB file name : ") #mol=read_pdb(f_name_pdb) mol=read_pdb('1chc.pdb') # file pointer for output fp_for_out=open('ramachandran_map.txt','w') for atom in mol: index=0 atom.gaps() phi, psi, omega=atom.torsions()[:3] for residue in atom.elements: fp_for_out.write('%10.2f %10.2f\n' %(phi[index],psi[index])) index+=1 fp_for_out.close()
1chc.pdb ATOM 1 N MET 1 66.104 56.583-35.505 1.00 0.00 1CHC 85 Calpha
Calpha-Calpha distance map (Python) (1/4) # homework 5, author: Shin, Seung-Woo # making distance matrix for Ca import math,string,os # Read PDB file input_file_name=raw_input("input the file name: ") fp_for_in=open(input_file_name,'r') fp_for_tmp=open("tmp.txt",'w') atom_count=0 # 2Pass processing because we don't know the number of Ca in advance # Step 1. Parsing Ca Atom part while 1: line=fp_for_in.readline() if not line: break else: token=string.split(line) if(token[0]=="atom" and token[2]=="ca"): line=token[0]+" "+token[5]+" "+token[6]+" "+token[7]+"\n" fp_for_tmp.write(line) atom_count=atom_count+1 fp_for_in.close() fp_for_tmp.close()
Calpha-Calpha distance map (Python) (2/4) # Step 2. read Ca X,Y,Z coordinates in array a_mat=[[0]*3 for i in xrange(atom_count)] fp_for_tmp=open("tmp.txt",'r') for i in range(0,atom_count): line=fp_for_tmp.readline() if not line: break else: token=string.split(line) a_mat[i][0]=float(token[1]) a_mat[i][1]=float(token[2]) a_mat[i][2]=float(token[3]) fp_for_tmp.close()
Calpha-Calpha distance map (Python) (3/4) # Define the matrix # Fill the distance matrix dist_matrix=[[0]*atom_count for i in xrange(atom_count)] for i in range(0,atom_count): for j in range(i,atom_count): dist=(a_mat[i][0]-a_mat[j][0])**2+(a_mat[i][1]-a_mat[j][1])**2+(a_mat[i][2]-a_mat[j][2])**2 dist=round(math.sqrt(dist)) if(dist>9): dist='.' dist_matrix[i][j]=dist dist_matrix[j][i]=dist else: dist_matrix[i][j]=int(dist) dist_matrix[j][i]=int(dist)
Calpha-Calpha distance map (Python) (4/4) fp_for_out=open("distance_mat.txt",'w') #=========================== #print distance matrix #=========================== # Write first row fp_for_out.write(" ") for i in range(1,atom_count+1): if i%10==0: fp_for_out.write("a") else: str0=i%10 str1=str(str0) fp_for_out.write(str1) fp_for_out.write("\n") #write dist_matrix for i in range(0,atom_count): i2=i+1 if i2<10: str1=" "+str(i2)+" " else: str1=str(i2)+" " fp_for_out.write(str1) for j in range(0,atom_count): str2=str(dist_matrix[i][j]) fp_for_out.write(str2) fp_for_out.write("\n") fp_for_out.close() # Delete temp file os.unlink("tmp.txt")
Fashion changes, style remains. - Gabrielle "Coco" Chanel Fashion is a form of ugliness so intolerable that we have to alter it every six months. - Oscar Wilde Every generation laughs at the old fashions, but follows religiously the new. - Henry David Thoreau
References http://www.python.org http://turing.cafe24.com/ http://www.dickbaldwin.com/tocpyth.htm http://www.diveintopython.org/toc/index.html http://pages.cs.wisc.edu/~fischer/cs538.s07/python.tutorial.pdf http://numpy.sourceforge.net/numdoc/html/numdoc.htm http://www.pasteur.fr/recherche/unites/sis/formation/python/index.html http://numpy.scipy.org/numpydoc/numpy.html http://numpy.scipy.org/numpy.pdf http://www.wellho.net/course/python.html http://www.pentangle.net/python/handbook/handbook.pdf http://docs.python.org/lib/lib.html http://rgruet.free.fr/pqr2.3.html http://www.awaretek.com/tutorials.html#begin http://www.biopython.org