Applications of Molecular Dynamics for Structural Analysis of Proteins

For a l l of these applications areas, extensive access to computers ... tion explores the effects of long range cutoffs on protein dynamics. This has...
0 downloads 0 Views 2MB Size
Chapter 8

Applications of Molecular Dynamics for Structural Analysis of Proteins and Peptides Bernard R. Brooks

Downloaded by EAST CAROLINA UNIV on January 3, 2018 | http://pubs.acs.org Publication Date: October 22, 1987 | doi: 10.1021/bk-1987-0353.ch008

Division of Computer Research and Technology, National Institutes of Health, Bethesda, MD 20892 A software system, GEMM (Generate, Emulate, and Manipulate Macromolecules), has been developed for the Star Technologies ST-100. This system allows the rapid computation of the energy and forces for use in molecular dynamics, minimization, or interactive modelling. A high degree of efficiency is obtained by the extensive development of microcode designed to run in a pipelined manner, which avoids the conventional problems of vectorization of the pair-interaction problem with a cutoff distance. The result of this implementation is that molecular dynamics simulations on the ST-100 are at least an order of magnitude more cost effective than on any other hardware system currently available. Several applications involving the use of this system are presented. The first application explores the effects of long range cutoffs on protein dynamics. This has been accomplished by carrying out 34 molecular dynamics simulations on carboxy-myoglobin. Each simulation is 150 picoseconds in length. For these calculations, both the distances at which the long range potentials are truncated and the method of truncation have been explored in a systematic manner. The second application demonstrates the utility of molecular dynamics for conformational exploration. This application involves the 3-dimensional structure determination of a glycoprotein N-terminal octapeptide T-cell epitope found on human ductal carcinoma (breast) cells. Over t h e l a s t s e v e r a l y e a r s , m o l e c u l a r dynamics has been becoming more u s e f u l a s a p r a c t i c a l t o o l f o r a wide v a r i e t y o f a p p l i c a t i o n s . These i n c l u d e e x p l o r i n g t h e temporal b e h a v i o r o f m a c r o m o l e c u l a r systems (1,2.), e x p l o r i n g energy decay ( 3 ) , s o l v e n t a n a l y s i s ( 4 ) , determining 3-dimensional s t r u c t u r e s u s i n g 2-dimensional n u c l e a r o v e r h a u s e r e f f e c t (NOE) d a t a from n u c l e a r magnetic r e s o n a n c e (NMR) experiments ( 5 ) , determining r e l a t i v e f r e e e n e r g i e s of s i m i l a r

This chapter not subject to U.S. copyright Published 1987 American Chemical Society

Jensen and Truhlar; Supercomputer Research in Chemistry and Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1987.

SUPERCOMPUTER

Downloaded by EAST CAROLINA UNIV on January 3, 2018 | http://pubs.acs.org Publication Date: October 22, 1987 | doi: 10.1021/bk-1987-0353.ch008

124

RESEARCH

systems ( 6 - 8 ) , c o n f o r m a t i o n a l sampling f o r s t r u c t u r e d e t e r m i n a t i o n (Brooks, B. R.; P a s t o r , R. W.; Carson, F. W. P r o c . N a t l . Acad. S c i . 1987, i n p r e s s . ) (9-11), and f o r c r y s t a l l o g r a p h i c R f a c t o r r e f i n e ment ( 1 2 ) . F o r a l l o f t h e s e a p p l i c a t i o n s a r e a s , e x t e n s i v e a c c e s s to computers w h i c h c a n r a p i d l y p e r f o r m m o l e c u l a r s i m u l a t i o n s i s essential. I n almost a l l p u b l i s h e d work i n t h i s f i e l d , t h e l e n g t h and q u a l i t y o f s i m u l a t i o n s i s b a l a n c e d by t h e a v a i l a b i l i t y o f computer r e s o u r c e s . F u r t h e r m o r e , as t h e s e s i m u l a t i o n methods become more r e f i n e d and r e l i a b l e , t h e t o t a l demand f o r s i m u l a t i o n c a p a b i l i t y w i l l i n c r e a s e many f o l d i n t h e next few y e a r s . One s o l u t i o n t o t h i s problem c u r r e n t l y b e i n g employed i s t h e s e t t i n g up o f l a r g e supercomputer f a c i l i t i e s from w h i c h many d i f f e r e n t r e s e a r c h e r s c a n benefit. A n o t h e r s o l u t i o n i s t o b u i l d s p e c i a l p u r p o s e hardware which i s s p e c i f i c a l l y d e s i g n e d t o p e r f o r m t h e time dominant s t e p ( t h e r e s t r i c t e d p a i r i n t e r a c t i o n problem) i n a m o l e c u l a r dynamics s i m u l a t i o n t o p r o v i d e s i m u l a t i o n c a p a b i l i t y i n a v e r y f a s t and c o s t e f f e c t i v e manner as i s b e i n g done f o r t h e FASTRUN p r o j e c t ( 1 3 ) . A t h i r d s o l u t i o n i s t o d e t e r m i n e t h e b e s t e x i s t i n g machine a r c h i t e c t u r e f o r t h e r e s t r i c t e d p a i r i n t e r a c t i o n and t o implement t h e simu l a t i o n on t h i s machine i n a h i g h l y o p t i m i z e d b u t g e n e r a l manner. An example o f t h i s t h i r d s o l u t i o n i s p r e s e n t e d i n t h i s paper which shows how i t i s p o s s i b l e t o a c h i e v e supercomputer speeds from a low c o s t l a b o r a t o r y computer. By p l a c i n g t h e t a s k on a l o c a l l a b comp u t e r , i t i s now a l s o p o s s i b l e t o d e v e l o p r e a s o n a b l e i n t e r a c t i v e m o l e c u l a r m o d e l l i n g t o o l s w h i c h u t i l i z e e n e r g i e s and f o r c e s i n r e a l time. The GEMM (Generate, Emulate, and M a n i p u l a t e M a c r o m o l e c u l e s ) system i s b e i n g used f o r almost a l l o f t h e s i m u l a t i o n work b e i n g p e r f o r m e d i n the M o l e c u l a r G r a p h i c s L a b o r a t o r y a t t h e N a t i o n a l I n s t i t u t e s o f H e a l t h (NIH). I n t h i s paper, p r e l i m i n a r y r e s u l t s o f two a p p l i c a t i o n s o f t h i s system a r e p r e s e n t e d , b o t h o f w h i c h would have been d i f f i c u l t t o c a r r y o u t on a s l o w e r machine. F u l l r e s u l t s o f t h e s e p r o j e c t s w i l l be p u b l i s h e d e l s e w h e r e . The f i r s t a p p l i c a t i o n e x p l o r e s t h e e f f e c t s o f l o n g range c u t o f f s on p r o t e i n dynamics. T h i s h a s been a c c o m p l i s h e d by c a r r y i n g o u t 34 m o l e c u l a r dynamics s i m u l a t i o n s on c a r b o x y - m y o g l o b i n . Each s i m u l a t i o n i s 150 p i c o seconds i n l e n g t h . F o r t h e s e c a l c u l a t i o n s , b o t h the d i s t a n c e s a t which t h e l o n g range p o t e n t i a l s a r e t r u n c a t e d and the method o f t r u n c a t i o n have been e x p l o r e d i n a s y s t e m a t i c manner. The second a p p l i c a t i o n d e m o n s t r a t e s the u t i l i t y o f m o l e c u l a r dynamics f o r c o n formational exploration. This a p p l i c a t i o n involves the structure determination of a g l y c o p r o t e i n N-terminal octapeptide T - c e l l e p i t o p e found on human d u c t a l c a r c i n o m a ( b r e a s t ) c e l l s . The

Star Technologies

ST-100

I n n o v a t i v e computer hardware d e s i g n s encourage new s o f t w a r e a l g o r i t h m s t h a t c a n p e r f o r m i n a v e r y e f f e c t i v e manner. One such example o f t h i s i s t h e development o f the m o l e c u l a r s i m u l a t i o n s o f t w a r e GEMM w h i c h has been d e v e l o p e d f o r t h e S t a r T e c h n o l o g i e s ST-100. GEMM i s a g e n e r a l p u r p o s e program t h a t c a n be a p p l i e d t o any t y p e o f m o l e c u l a r system. A t b e s t t h e ST-100 i s c a p a b l e o f 100 m i l l i o n f l o a t i n g p o i n t o p e r a t i o n s p e r second (MFLOPS); however f o r r e a s o n a b l y l a r g e macro-

Jensen and Truhlar; Supercomputer Research in Chemistry and Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1987.

Downloaded by EAST CAROLINA UNIV on January 3, 2018 | http://pubs.acs.org Publication Date: October 22, 1987 | doi: 10.1021/bk-1987-0353.ch008

8.

BROOKS

Structural

Analysis

of Proteins and Peptides

125

m o l e c u l a r systems, t h e p e r f o r m a n c e o f t h i s s o f t w a r e runs a t about 90% o f t h e speed o f an o p t i m i z e d Cray X-MP v e r s i o n o f s i m i l a r s o f t ware, and even f a s t e r than t h e Cray X-MP f o r s m a l l e r systems where vectorization i s d i f f i c u l t . The throughput e f f i c i e n c y o f t h e s o f t ware on t h e ST-100 i s r o u g h l y 50% ( i . e . 50 MFLOPS). T h i s h i g h d e g r e e o f e f f i c i e n c y i s o b t a i n e d by t h e e x t e n s i v e development o f m i c r o c o d e d e s i g n e d t o r u n i n a p i p e l i n e d manner, w h i c h a v o i d s t h e c o n v e n t i o n a l problems o f v e c t o r i z a t i o n o f t h e p a i r - i n t e r a c t i o n problem w i t h a c u t o f f d i s t a n c e . The r e s u l t o f t h i s implementation i s t h a t m o l e c u l a r dynamics s i m u l a t i o n s on t h e ST-100 a r e a t l e a s t an o r d e r o f magnitude more c o s t e f f e c t i v e than on any o t h e r h a r d ware system c u r r e n t l y a v a i l a b l e . The GEMM s o f t w a r e on t h e ST-100 i s n o t a s t a n d - a l o n e package, and i t r e q u i r e s a f r o n t - e n d s i m u l a t i o n s o f t w a r e package t h a t r u n s on the h o s t t o p r o v i d e d a t a and t o send command r e q u e s t s . I t was d e s i g n e d and w r i t t e n w i t h CHARMM ( C h e m i s t r y a t HARvard Macrom o l e c u l a r M e c h a n i c s ) (14) a s t h e p r i m a r y f r o n t - e n d , b u t a d d i t i o n a l s o f t w a r e packages, such a s AMBER ( 1 5 ) , have s u b s e q u e n t l y been m o d i f i e d t o d r i v e GEMM. The ST-100 i s n o t a t y p i c a l a r r a y p r o c e s s o r . F i g u r e 1 shows the d a t a f l o w t h r o u g h t h e f u n c t i o n a l u n i t s o f the machine f o r a routine application. T h e r e a r e seven main f u n c t i o n a l u n i t s ; (1) The C o n t r o l P r o c e s s o r r u n s a FORTRAN-like language w h i c h c o n t r o l s a l l o f t h e f u n c t i o n a l elements w i t h i n t h e ST-100. (2) The i n p u t / o u t p u t (I/O) subsystem h a n d l e s a l l d a t a t r a n s f e r between t h e h o s t and t h e ST-100 and i s t r a n s p a r e n t t o the u s e r . (3) The main memory i s t y p i c a l l y 8 megabytes and i s b y t e a d d r e s s a b l e and has a d a t a t r a n s f e r r a t e o f 800 m e g a b i t s / s e c o n d . (4) The S t o r a g e Move P r o c e s s o r (SMP) moves d a t a between the main and cache memories, c o n v e r t s d a t a between h o s t and i n t e r n a l f l o a t i n g p o i n t f o r m a t , and p e r f o r m s i n t e g e r and B o o l e a n o p e r a t i o n s ( b i t p a c k i n g , e t c . ) . (5) The cache memory c o n s i s t s o f t h r e e banks o f memory e a c h 16384 3 2 - b i t words l o n g . Each cache bank c a n be s p l i t i n t o h a l v e s and e i t h e r h a l f c a n be a s s i g n e d t o t h e SMP o r t o the A r i t h m e t i c C o n t r o l P r o c e s s o r (ACP). T h i s allows the double b u f f e r i n g o f data t r a n s f e r s between main memory and t h e cache memory. Data t r a n s f e r s between the c o m p u t a t i o n a l elements and d i f f e r e n t cache banks c a n o c c u r a t the same time. (6) The ACP p r o c e s s e s the m i c r o c o d e t h a t d r i v e s the c o m p u t a t i o n a l e l e m e n t s . I t a l s o performs i n t e g e r a r i t h m e t i c , cache bank a d d r e s s i n g , and c o n d i t i o n a l b r a n c h i n g . (7) The comput a t i o n a l elements a r e i n d e p e n d e n t and communicate t h r o u g h a d a t a interchange device. The add and m u l t i p l y elements have t h r e e s t a g e s so t h a t a new o p e r a t i o n may be s t a r t e d on e v e r y c y c l e , b u t i t r e q u i r e s three c y c l e s f o r a s i n g l e o p e r a t i o n t o complete. Since the machine c y c l e i s 40 nanoseconds, each element c a n p e r f o r m a t 25 MFLOPS. The d i v i d e / s q u a r e - r o o t element i s n o t p i p e l i n e d and 15 cycles are required f o r a d i v i d e . When an a p p l i c a t i o n r e q u i r e s a d i v i d e , any o t h e r a r i t h m e t i c o p e r a t i o n c a n p r o c e e d c o n c u r r e n t l y i n the o t h e r c o m p u t a t i o n a l e l e m e n t s . By u s i n g t h e d a t a i n t e r c h a n g e i n a v a r i e t y o f ways i t i s p o s s i b l e t o program t h i s machine i n t h e manner o f a c o n v e n t i o n a l a r r a y p r o c e s s o r o r a s a p i p e l i n e p r o c e s s o r w i t h up t o about 15 c o n c u r r e n t c a l c u l a t i o n s . T h i s i s made p o s s i b l e because t h e r e a r e no v e c t o r r e g i s t e r s a s t h e r e a r e on most v e c t o r machines, and a l l memory a c c e s s c a n be c o n t r o l l e d by the u s e r .

Jensen and Truhlar; Supercomputer Research in Chemistry and Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1987.

Jensen and Truhlar; Supercomputer Research in Chemistry and Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1987.

si n

W

H M

C

η ο S

M

c/5 G

fcO

Downloaded by EAST CAROLINA UNIV on January 3, 2018 | http://pubs.acs.org Publication Date: October 22, 1987 | doi: 10.1021/bk-1987-0353.ch008

8.

BROOKS

Structural

Analysis

of Proteins and Peptides

127

Goals. The p r i m a r y g o a l i n d e v e l o p i n g GEMM was t o p r o v i d e a v e r y f a s t g e n e r a l purpose m o l e c u l a r s i m u l a t i o n package. T h i s r e q u i r e s the r a p i d c o m p u t a t i o n o f e n e r g i e s and f o r c e s f o r t h e u s e o f m o l e c ­ u l a r dynamics, energy m i n i m i z a t i o n , and i n t e r a c t i v e m o d e l l i n g . I t was d e c i d e d e a r l y i n development t h a t t h e r e would be no s a c r i f i c e o f t h e b e s t f u n c t i o n a l forms t h a t have been used i n CHARMM and o t h e r m o l e c u l a r s i m u l a t i o n programs (14, 1 5 ) . F o r example, a g r e a t d e a l o f s o f t w a r e development e f f o r t c o u l d have been a v o i d e d i f s i m p l e t r u n c a t i o n o f l o n g range f o r c e s ( e l e c t r o s t a t i c and v a n d e r Waals) had been s u f f i c i e n t . I n a d d i t i o n , such a r e s t r i c t i o n would a l l o w t h e program t o r u n s i g n i f i c a n t l y f a s t e r (Hagstrom, R., Argonne N a t i o n a l L a b o r a t o r i e s , p e r s o n a l communication, 1986.). Nevertheless, the p r e s e n c e o f l a r g e d i s c o n t i n u i t i e s i n t h e f o r c e s makes m i n i m i z a ­ t i o n s v e r y d i f f i c u l t , and h a s s e v e r a l d e l e t e r i o u s e f f e c t s on m o l e c ­ u l a r dynamic s i m u l a t i o n s ( L e v i t t , Μ . , Weizmann I n s t i t u t e , p e r s o n a l communication, 1986.). As t h e program i s c u r r e n t l y w r i t t e n , t h e r e a r e many ways i n w h i c h t h e l o n g range f o r c e terms c a n be t e r m i n a t e d as a f u n c t i o n o f d i s t a n c e . One p e r c e i v e d problem w i t h t h e ST-100 f o r m o l e c u l a r s i m u l a t i o n s i s t h a t a l l a r i t h m e t i c f l o a t i n g p o i n t o p e r a t i o n s a r e done i n s i n g l e p r e c i s i o n (24- o r 3 2 - b i t m a n t i s s a , depending on o p e r a t i o n ) , b u t some o p e r a t i o n s r e q u i r e higher p r e c i s i o n , p a r t i c u l a r l y those i n v o l v e d w i t h a c c u m u l a t i n g f o r c e components. T h i s problem has been overcome by t h e c r e a t i o n o f s i x d o u b l e p r e c i s i o n macros ( 4 8 - b i t m a n t i s s a ) w h i c h a r e employed wherever h i g h e r p r e c i s i o n i s n e c e s s a r y . These macros u s e s e v e r a l s i n g l e p r e c i s i o n o p e r a t i o n s t o compute a s i n g l e high precision r e s u l t . As a r e s u l t , a l l i n d i v i d u a l energy and f o r c e components a r e computed i n low p r e c i s i o n , b u t a l l a c c u m u l a t i o n s a r e done i n h i g h e r p r e c i s i o n . T h i s g i v e s most o f the b e n e f i t s o f h i g h p r e s i c i o n , w i t h o u t t h e extreme d i f f i c u l t i e s t h a t would e n t a i l i n t h e development o f a f u l l d o u b l e p r e c i s i o n v e r s i o n . A n o t h e r g o a l o f t h e development o f GEMM was t o a l l o w v e r y l a r g e systems t o be s t u d i e d . T h i s was a c c o m p l i s h e d by u s i n g a compact form o f the nonbond l i s t where, on t h e a v e r a g e , each element r e ­ q u i r e s o n l y one b y t e . T h i s c o n t r a s t s t o t h e two b y t e s employed by CHARMM f o r a VAX and t h e e i g h t b y t e s t h a t a r e employed on a C r a y . Thus f o r a t y p i c a l ST-100 i n s t a l l a t i o n ( w i t h 2 megawords o f memory), a system w i t h 8000 atoms and up t o 7 m i l l i o n nonbond p a i r s c a n be s i m u l a t e d (4000 atoms w i t h h i g h p r e c i s i o n ) . W i t h t h e 4x-cache memory o p t i o n , systems w i t h 16000 atoms c a n be s t u d i e d i n h i g h precision. The l i m i t o f 7 m i l l i o n nonbond p a i r s a l l o w s v e r y dense p e r i o d i c systems t o be s t u d i e d such a s a p r o t e i n i n a box o f water, or a l i p i d b i l a y e r . A n o t h e r programming c o n s i d e r a t i o n i s t h e e f f e c t i v e u s e o f t h e I/O c h a n n e l between t h e ST-100 and t h e h o s t . W i t h a VAX h o s t , t h e d a t a t r a n s f e r r a t e s a r e l e s s than 1 m e g a b i t / s e c o n d . This precludes the p o s s i b i l i t y o f h a v i n g a f a s t system where t h e t h r e e - d i m e n s i o n a l c o o r d i n a t e s a r e t r a n s f e r r e d from the h o s t t o t h e ST-100 and t h e sub­ s e q u e n t l y computed e n e r g i e s and f o r c e s r e t r i e v e d from t h e ST-100 on every step. The GEMM package has been d e s i g n e d t o m i n i m i z e synchronous d a t a t r a n s f e r s . I n an i d e a l s i t u a t i o n , t h e ST-100 runs w i t h o u t d i r e c t i o n from t h e h o s t , and d a t a t h a t need t o be saved a r e streamed t o the h o s t w h i l e t h e ST-100 c o n t i n u e s t o p r o d u c e new d a t a .

Jensen and Truhlar; Supercomputer Research in Chemistry and Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1987.

128

SUPERCOMPUTER RESEARCH

Downloaded by EAST CAROLINA UNIV on January 3, 2018 | http://pubs.acs.org Publication Date: October 22, 1987 | doi: 10.1021/bk-1987-0353.ch008

Features. GEMM i s w r i t t e n i n a h o s t - i n d e p e n d e n t manner and i t h a s been r u n w i t h an A p o l l o , a VAX, and a MicroVAX I I as a h o s t . GEMM can c u r r e n t l y p e r f o r m t h e f o l l o w i n g o p e r a t i o n s : perform molecular dynamics, p e r f o r m energy m i n i m i z a t i o n s , compute t h e energy and f o r c e s f o r a s t r u c t u r e , and u p d a t e t h e nonbond l i s t (nonbond l i s t s are u s u a l l y automatic f o r the other o p e r a t i o n s ) . In addition, a wide v a r i e t y o f I/O sequences a r e p o s s i b l e , s u c h a s what i s needed f o r i n t e r a c t i v e m o d e l l i n g work. For l a r g e systems, t h e c o m p u t a t i o n o f t h e energy c a n r e q u i r e 98% o r more o f t h e t o t a l c o m p u t a t i o n a l e f f o r t , and w i t h i n t h e energy c o m p u t a t i o n , the p a i r i n t e r a c t i o n energy (van d e r Waals and e l e c t r o s t a t i c terms) c a n r e p r e s e n t more t h a n 90% o f t h e t o t a l effort. F o r t h i s r e a s o n , t h e a s p e c t s o f t h e program t h a t d e a l w i t h the g e n e r a t i o n o f t h e nonbond l i s t and t h e c o m p u t a t i o n o f the n o n bond energy a r e o f c r i t i c a l i m p o r t a n c e . The nonbond l i s t g e n e r a t i o n p a i r s e a r c h i n g s e c t i o n f i n d s a l l p a i r s o f atoms w i t h i n a g i v e n c u t o f f d i s t a n c e , u s u a l l y 8 angstroms o r more. The i n n e r l o o p o f t h i s code h a s 5 c y c l e s (0.2 m i c r o s e c o n d s / p a i r ) and p e r f o r m s w i t h 80% e f f i c i e n c y (80 MFLOPS). I t a l l o w s atoms t o be f i x e d ( f r o z e n ) , and t h e r e i s a l s o a f a c i l i t y t o s e l e c t by g r o u p s , so t h a t a l l atom p a i r s between two groups a r e e i t h e r i n c l u d e d o r e x c l u d e d a s a whole. T h i s code works i n a d o u b l e b u f f e r e d mode ( a s y n c h r o n o u s l y ) w i t h t h e l i s t merging o p e r a t i o n . The l i s t merging o p e r a t i o n removes from t h e c l o s e c o n t a c t l i s t (from t h e p a i r s e a r c h i n g s t e p ) t h e nonbond e x c l u s i o n s , w h i c h i n c l u d e s b o n d i n g (1-2) and a n g l e (1-3) i n t e r a c t i o n s , t o g e n e r a t e t h e f i n a l nonbond pair l i s t . The f i n a l l i s t i s packed i n a b y t e f o r m a t . The l i s t merging macro h a s 8 c y c l e s i n t h e i n n e r l o o p , b u t o n l y a c t s on t h o s e elements a c c e p t e d by the p a i r s e a r c h i n g code, so t h e c o s t o f t h i s o p e r a t i o n i s c o m p l e t e l y h i d d e n e x c e p t i n c a s e s where the n o n bond c u t o f f d i s t a n c e i s l a r g e i n r e l a t i o n t o t h e s i z e o f the m o l e c u l a r system. T h i s code i s e x t r e m e l y messy and i t u s e s unannounced f e a t u r e s o f t h e ST-100 t o a c c o m p l i s h l i s t p r o c e s s i n g . There i s a l s o s p e c i a l t r e a t m e n t p r o v i d e d f o r 1-4 i n t e r a c t i o n s so t h a t s e p a r a t e van d e r Waals p a r a m e t e r s c a n be s p e c i f i e d and t h e e l e c t r o s t a t i c c o n t r i b u t i o n s c a n be r e d u c e d . The nonbond energy e v a l u a t i o n c o n s i s t s o f t h r e e macros. Since e a c h i n d i v i d u a l macro r e q u i r e s s e v e r a l thousand machine c y c l e s t o complete, t h e r e i s no expense i n b r e a k i n g t h e problem up i n t o s e v e r a l p a r t s b e c a u s e t h e l o a d i n g o f macros i s a s y n c h r o n o u s w i t h t h e i r execution. The f i r s t macro computes t h e r e c i p r o c a l o f t h e s q u a r e o f t h e d i s t a n c e between i n t e r a c t i n g atoms ( l / r * * 2 ) , removes atom p a i r s beyond t h e c u t o f f d i s t a n c e , and p r o c e s s e s t h e minimum image o p t i o n . Because t h e c o s t o f p r o c e s s i n g the minimum image c o n v e n t i o n i s h i d d e n b e h i n d t h e d i v i d e , t h i s code i s always p r e s e n t and f o r c a l c u l a t i o n s where t h e r e a r e no p e r i o d i c b o u n d a r i e s , a v e r y l a r g e box s i z e i s u s e d . The u s e o f m u l t i p l e macros s i m p l i f i e s t h e c o d i n g and a l l o w s r e p l a c e m e n t o f d i f f e r e n t s e c t i o n s t o a l l o w f o r d i f f e r e n t options. F o r example, t h e r e a r e two v e r s i o n s o f t h e m i d d l e macro, one t h a t a l l o w s a s w i t c h i n g f u n c t i o n f o r l o n g range i n t e r a c t i o n t e r m i n a t i o n , and one t h a t u s e s a " s h i f t e d " p o t e n t i a l . T h e r e a r e a l s o two v e r s i o n s o f t h e l a s t macro. One t h a t a c c u m u l a t e s f o r c e s i n s i n g l e p r e c i s i o n and one t h a t a c c u m u l a t e s f o r c e s i n d o u b l e precision. W i t h i n the m i d d l e macros, t h e r e i s code t h a t computes

Jensen and Truhlar; Supercomputer Research in Chemistry and Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1987.

Downloaded by EAST CAROLINA UNIV on January 3, 2018 | http://pubs.acs.org Publication Date: October 22, 1987 | doi: 10.1021/bk-1987-0353.ch008

8.

BROOKS

Structural

Analysis

of Proteins and

Peptides

129

the energy w i t h a c o n s t a n t d i e l e c t r i c o r w i t h a d i e l e c t r i c t h a t i s p r o p o r t i o n a l to the d i s t a n c e . Having many o p t i o n s w i t h no c o s t i n e f f i c i e n c y makes the program u s e f u l f o r a wide c l a s s o f p r o b l e m s , and f o r a w i d e r d i s t r i b u t i o n o f u s e r s . T h e r e a r e a l s o t h e i n t e r n a l energy terms, t h o s e t h a t a c c o u n t f o r t h e c o n n e c t i v i t y of groups o f atoms (bonds, a n g l e s , d i h e d r a l , and improper d i h e d r a l s ) . W i t h i n t h e s e terms, a l l f o r c e s can be accumulated i n h i g h p r e c i s i o n , and t h e r e i s d o u b l e b u f f e r i n g b o t h w i t h i n and between terms. F o r example, when t h e l a s t a n g l e energy b u f f e r i s b e i n g computed, t h e f i r s t d i h e d r a l b u f f e r i s b e i n g l o a d e d i n t o cache memory, thus c o m p l e t e l y h i d i n g the c o s t o f l o a d i n g the d i h e d r a l d a t a a r r a y s and t h e g a t h e r i n g o f t h e c o o r d i n a t e s . For s m a l l systems where t h e s e terms can dominate t h e t o t a l c o m p u t a t i o n a l e f f o r t , t h i s can be i m p o r t a n t , but u s u a l l y , t h e s e terms a r e n o t time i n t e n s i v e . N e v e r t h e l e s s , t h i s a c c o u n t s f o r much o f the code (and c o d i n g e f f o r t ) o f GEMM. (See T a b l e I.) Table I.

T i m i n g s f o r Each S e c t i o n o f Code

of Code S e c t i o n macros Nonbond l i s t s e a r c h i n g 1 Nonbond l i s t f o r images 1 Nonbond l i s t merging 1 Nonbond energy 3 Bond energy 3 A n g l e energy 9 D i h e d r a l energy 7 Improper t o r s i o n energy 9 Hydrogen bond energy 26 Harmonic c o n s t r a i n t s 2 SHAKE c o n s t r a i n t s 3

// o f c y c l e s singledoubleprecision precision 5 11 8 34 43 51 41 144 159 130 150 178 198 294 254 23 18 40 35

time/element (μ seconds) 0.2 0.4 0.3 1.4 1.6 5.8 5.2 7.1 10.2 0.7 1.4

Performance. The t o t a l time f o r 500 s t e p s o f dynamics and 25 nonbond u p d a t e s f o r t h e s t a n d a r d CHARMM benchmark (Brunger, A. T., H a r v a r d U n i v e r s i t y , p e r s o n a l communication, 1985.), a B-DNA e l e v e n mer d u p l e x w i t h 706 atoms and a 11.5 angstrom nonbond c u t o f f (77000 nonbond p a i r s ) i s found i n T a b l e I I . Table I I .

CHARMM Performance

on a V a r i e t y o f

Machine Time ( s e c ) A p o l l o DSP-160 d o u b l e p r e c i s i o n 22932 A p o l l o DSP-160 s i n g l e p r e c i s i o n 14904 VAX 11/780 w/FPA d o u b l e p r e c i s i o n 11650 IBM 3090 w/vector d o u b l e p r e c i s i o n 226 C r a y - l S w/vector 186 Cray-lS w/vector(*) 174 Cray X-MP/48 no g a t h e r / s c a t t e r ( * ) 140 ST-100 d o u b l e p r e c i s i o n ( * ) 89.9 ST-100 s i n g l e p r e c i s i o n ( * ) 76.8 Cray X-MP/48 w i t h g a t h e r / s c a t t e r ( * ) 70 (*) does n o t i n c l u d e s e t u p time

Machines

R a t i o t o ST-100 299 194 152 2.94 2.42 2.27 1.82 1.17 1.00 0.91

Jensen and Truhlar; Supercomputer Research in Chemistry and Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1987.

SUPERCOMPUTER RESEARCH

130

The t i m i n g s r e p r e s e n t t h o s e u s i n g the b e s t code f o r each p a r t i c u l a r machine. F o r t h e system chosen, the v e c t o r l e n g t h s a r e l o n g e r t h a n u s u a l on a v e r a g e . F o r s m a l l systems t h a t do not v e c t o r i z e w e l l t h e ST-100 does even b e t t e r i n r e l a t i o n t o the o t h e r v e c t o r machines because o f i t s p i p e l i n e c o d i n g and heavy u s e o f d o u b l e b u f f e r i n g a c r o s s t h e d i f f e r e n t energy terms. For p r o t e i n s w i t h a polar hydrogen r e p r e s e n t a t i o n and an 8 angstrom c u t o f f , s i n g l e p r e c i s i o n , and no e x t e r n a l water, t h e t i m i n g s a r e g i v e n i n T a b l e I I I .

Downloaded by EAST CAROLINA UNIV on January 3, 2018 | http://pubs.acs.org Publication Date: October 22, 1987 | doi: 10.1021/bk-1987-0353.ch008

Table I I I .

Timings P e r I t e r a t i o n f o r D i f f e r e n t

Energy o n l y System (seconds) D u c t a l carcinoma o c t a p e p t i d e 0.,014 BPTI (580 atoms) 0.,070 Lysozyme (1264 atoms) 0,,162 Carboxy-myoglobin (1532 atoms) 0.,174 Hemoglobin (5478 atoms) 0.,518

Proteins

With l i s t updates every 20th s t e p (seconds) N/A 0.072 0.171 0.189 0.674

When one c o n s i d e r s t h a t t h e ST-100 i s a $275k machine and the h a l f speed ST-50 has a l i s t p r i c e o f $125k, t h e r e a l p o t e n t i a l f o r t h e s e machines becomes a p p a r e n t . Future D i r e c t i o n s . The GEMM package i s s t i l l under development. T h e r e a r e s e v e r a l a r e a s where new c a p a b i l i t i e s a r e f o r t h c o m i n g . One f e a t u r e j u s t completed i s t h e i n c o r p o r a t i o n o f p e r i o d i c boundary c o n d i t i o n s so t h a t systems i n a water box may be s t u d i e d . Another f e a t u r e t h a t h a s j u s t been completed and i s now b e i n g t e s t e d i s t h e a b i l i t y t o p e r f o r m L a n g e v i n dynamics, so t h a t systems c a n be c o u p l e d to a h e a t b a t h ( 1 6 ) . An a r e a o f development i s t h e i n c l u s i o n o f code n e c e s s a r y f o r f r e e energy p e r t u r b a t i o n c a l c u l a t i o n s f o r i n t e r n a l , v a n d e r Waals and e l e c t r o s t a t i c energy terms ( 6 ) . A n o t h e r a r e a o f development i s the o p t i o n t o u s e lookup t a b l e s f o r b o t h nonbond e n e r g i e s and f o r c o n s t r a i n t energy terms s u c h a s a s o l v e n t boundary energy term ( 1 6 ) . O t h e r f e a t u r e s and o p t i o n s a r e r o u t i n e l y added a s needed t o a i d i n the e x e c u t i o n o f c u r r e n t p r o j e c t s underway w i t h i n the NIH. A n o t h e r a r e a o f f u t u r e p r o m i s e i s t h e i n t e r f a c e o f GEMM w i t h the FASTRUN d e v i c e c u r r e n t l y b e i n g d e v e l o p e d a t Columbia U n i v e r s i t y under t h e d i r e c t i o n o f P r o f e s s o r Cyrus L e v i n t h a l ( 1 3 ) . The FASTRUN d e v i c e i s a c o m p u t a t i o n a l u n i t which i s d i r e c t l y i n t e r f a c e d t o a ST-100 t h r o u g h a Dataram memory u n i t and i s e x p e c t e d t o p e r f o r m a t 700 MFLOPS. I t s h o u l d i n c r e a s e t h e throughput o f GEMM by a f a c t o r of 5 t o 10 o v e r i t s c u r r e n t p e r f o r m a n c e . Consideration for i n t e r f a c i n g w i t h FASTRUN h a s been k e p t i n t h e c u r r e n t d e s i g n o f GEMM, and once FASTRUN i s completed, a w o r k i n g v e r s i o n s h o u l d s u r f a c e w i t h i n t h r e e months. Conclusions. GEMM i s a program f o r m o l e c u l a r s i m u l a t i o n s w r i t t e n f o r t h e ST-100 ( o r ST-50) w h i c h p e r f o r m s a l m o s t a t Cray X-mp speeds f o r l a r g e systems. As such, GEMM r u n n i n g on t h e ST-100 ( o r ST-50) i s an o r d e r o f magnitude c o s t e f f e c t i v e i n comparison w i t h any o t h e r hardware system c u r r e n t l y a v a i l a b l e . GEMM i s w r i t t e n i n a p i p e l i n e manner t h e r e b y a v o i d i n g t h e problems o f v e c t o r i z a t i o n found on a

Jensen and Truhlar; Supercomputer Research in Chemistry and Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1987.

Downloaded by EAST CAROLINA UNIV on January 3, 2018 | http://pubs.acs.org Publication Date: October 22, 1987 | doi: 10.1021/bk-1987-0353.ch008

8.

BROOKS

Structural

Analysis

of Proteins and Peptides

131

c o n v e n t i o n a l machine. Most o f the program i s w r i t t e n i n m i c r o c o d e (^90%) and a s such i t i s f a s t , b u t r a t h e r i n f l e x i b l e . Itis diffi­ c u l t t o change t h e f u n c t i o n a l forms o r t o add new ones, and major new f e a t u r e s r e q u i r e c a r e f u l development. F u r t h e r m o r e , t h e program d e s i g n i s n o t a p p r o p r i a t e f o r any o t h e r machine o t h e r than t h e ST-100 o r ST-50. As such, t h e ST-100 i s n o t t h e b e s t machine f o r the development o f new methods, b u t i t i s b e t t e r used as a t o o l i n a p p l i c a t i o n a r e a s where the methods a r e p r o v e n . An advantage o f t h e ST-100 i s t h a t i t i s a l a b o r a t o r y machine. The problems a s s o c i a t e d w i t h r u n n i n g s i m u l a t i o n s on a remote s i t e do not a p p l y . As a l a b o r a t o r y t o o l , t h e ST-100 i s i d e a l f o r i n t e r a c ­ t i v e u s e such as m o l e c u l a r m o d e l l i n g w i t h r e a l time energy d e t e r m i ­ n a t i o n o r energy m i n i m i z a t i o n s . The problems a s s o c i a t e d w i t h low p r e c i s i o n a t t h e hardware l e v e l have been overcome, and GEMM has been d e s i g n e d t o t r e a t l a r g e systems w i t h up t o 1600 atoms i n h i g h p r e c i s i o n and up t o 7 m i l l i o n nonbonding atom p a i r s f o r a machine w i t h two megawords o f main memory. The c u r r e n t i m p l e m e n t a t i o n i s b e s t w i t h CHARMM as a f r o n t - e n d , though o t h e r s o f t w a r e packages have been m o d i f i e d t o d r i v e GEMM. The o n l y h o s t machines t e s t e d have been a VAX, a MicroVAX I I , and the A p o l l o w o r k s t a t i o n s . GEMM i s a v a i l a b l e f o r d i s t r i b u t i o n . The

E f f e c t s o f Long Range C u t o f f s on P r o t e i n Dynamics

Throughout t h e f i e l d o f m o l e c u l a r s i m u l a t i o n s t h e r e a r e many ways i n which t h e t r u n c a t i o n o f l o n g range f o r c e s , p a r t i c u l a r l y e l e c t r o ­ s t a t i c forces, i s achieved. T h e r e have been s e v e r a l p a p e r s d e t a i l ­ i n g t h e e f f e c t s o f l o n g range t r u n c a t i o n f o r p o l a r and i o n i c s o l ­ v e n t s u s i n g a n a l y t i c and s t a t i s t i c a l methods ( 1 7 ) . T h e r e have a l s o been examples o f c a s e s where problems found i n t h e a n a l y s i s o f s i m u l a t i o n s have r e s u l t e d i n a b e t t e r t r e a t m e n t o f t h e l o n g range p o t e n t i a l ( L e v i t t , Μ . , t o be p u b l i s h e d ) , b u t most o f t h e s e a r e imbedded i n p a p e r s n o t d e a l i n g s p e c i f i c a l l y w i t h l o n g range t r u n c a ­ tion. T h e r e has a l s o been t h e tendency t o d i m i n i s h some a d v e r s e a f f e c t s , s u c h a s d i s c o n t i n u i t i e s i n t h e f o r c e , by c o u p l i n g t h e e n t i r e system t o a h e a t b a t h (18) r a t h e r t h a n remove t h e d i s c o n t i n u ­ ity. I n a l l o f t h e s e c a s e s , t h e r e a r e t r a d e - o f f s between e f f i c i e n c y , a c c u r a c y , ease o f programming, v e c t o r i z a b i l i t y , and j u s t i f i a b i l i t y . Due t o t h e c o m p l e x i t y o f p r o t e i n - p r o t e i n i n t e r a c t i o n s , i t i s n o t f e a s i b l e t o a c c u r a t e l y p r e d i c t t h e e f f e c t s o f l o n g range t r u n c a t i o n by methods o t h e r t h a n d i r e c t s i m u l a t i o n . Methods. To e x p l o r e t h e e f f e c t s o f l o n g range t r u n c a t i o n on m o l e c ­ u l a r dynamics s i m u l a t i o n s o f p r o t e i n s , a s e r i e s o f s i m u l a t i o n s have been p e r f o r m e d . The two main q u e s t i o n s which have been a d d r e s s e d are: 1) What i s t h e b e s t f u n c t i o n a l form f o r t r u n c a t i n g l o n g range forces? 2) What i s t h e b e s t d i s t a n c e t o e f f e c t t r u n c a t i o n ? S e v e r a l a d d i t i o n a l q u e s t i o n s have been a d d r e s s e d s u c h a s : 3) What i s t h e e f f e c t o f a s h o r t s w i t c h i n g f u n c t i o n ? 4) What i s t h e r e l a t i o n s h i p between a c o n s t a n t d i e l e c t r i c and a d i s t a n c e dependant d i e l e c t r i c f o r pseudo-vacuum s i m u l a t i o n s ? 5) What i s t h e e f f e c t o f i n c r e a s i n g t h e d i e l e c t r i c c o n s t a n t ?

Jensen and Truhlar; Supercomputer Research in Chemistry and Chemical Engineering ACS Symposium Series; American Chemical Society: Washington, DC, 1987.

SUPERCOMPUTER RESEARCH

132

Four methods f o r t r u n c a t i o n have been e x p l o r e d . (a) s i m p l e atom-atom t r u n c a t i o n q q E

ii

=

F ( r )

Downloaded by EAST CAROLINA UNIV on January 3, 2018 | http://pubs.acs.org Publication Date: October 22, 1987 | doi: 10.1021/bk-1987-0353.ch008

J

These i n c l u d e :

F(r) =1

r < r

F(r) =0

r ^ r _ cut

(

—^27 » 4πεΓ^ ;

7

1

)

T h i s f u n c t i o n a l form i s e a s i e s t t o program and t o v e c t o r i z e . T h e r e a r e d i s c o n t i n u i t i e s i n t h e energy and f o r c e s which make t h i s i n a p ­ p r o p r i a t e f o r m i n i m i z a t i o n s , and causes problems f o r dynamics p a r ­ t i c u l a r l y i f solvent i s included. (b) an atom-atom s w i t c h i n g f u n c t i o n o v e r a range (14) = 1

r>

&fl t.