Documentation set for the MOIL package. Description of programs: BOAT.PRG (BOnd Angles and Torsions) ............................... 2 CCRD.PRG (Coordinate conversion) .................................. 4 CHAIN.PRG (Simulated annealing to find minimum energy path) ........ 5 CON.PRG (Forms molecules from a data base) ....................... 13 DYNA.PRG (Molecular dynamics program) ............................. 20 MINI.PRG (Minimization program) ................................... 25 RMS.PRG (Root Mean Square deviation between two structures) ...... 27 Description of data structures: CONN.FIL (Connectivity of a molecule) .............................. 29 MONO.FIL (Monomer connectivity file) ............................... 33 POLY.FIL (Sequence information used by CON.PRG) .................... 35 PROP.FIL (Property file) ........................................... 36 SUB.LIST (List of subroutines and their functions) ................. 39 TUTORIAL.SHORT (A good place to start) ............................. 43 WARNING IT IS PROBABLY SAFER TO LOOK FOR EXACT SYNTAX IN THE moil.tests DIRECTORY SINCE THE LAST IS UP to DATE WHILE THESE DOCUMENTATION ARE NOT. ALL INPUT FILE SHOULD WORK AS USUAL. HOWEVER SOME NEW FEATURES MAY BE MISSING. 1) Relatively important cosmetic changes includes the removal of the unit=[int] requirement in input. The program is doing its own file handling 2) As of now there are no documentation for the freee and umbrella programs see moil.tests for examples. In general the syntax is quite similar to dyna 3) The "mini" program was splitted into to mini_pwl (Powell conjugate gradient minimization) and mini_tn (truncated newton method) BOAT.PRG Here is how the current values of the BOnds Angles and Torsions (BOAT) can be calculated: Needed files are (i) connectivity (ii) coordinates. Out files include the standard output and another file with BOAT listing. The input below assumes that the connectivity is in VAL.WCON and the coordinates in CMIN.CRD. The coordinates type is defined by the line coor CHAR They can be in CHARMM format (only one set), DYNA format and PATH format. The last two should come with a range of required structures, e.g., coor PATH lpst=2 lpen=10 Thus, find all BOAT's starting from structure 2 to structure 10 in the PATH formatd file, provided by the "file rcrd ..." command. Output file is in VAL.BOAT (Note that the connectivity file is used directly, i.e. torsions that do not appear in the conn file will not be calculated. This may suggest tampering in PROP or WCON file to make sure that torsion you are interested in will be printted out and other that your are not will not. ~ BOAT TEST ~ ~debug file rcon name=(VAL.WCON) unit=10 read file rcrd name=(CMIN.CRD) unit=11 read file boat name=(VAL.BOAT) unit=13 wovr coor CHAR action action initiates appropriate(?) activity. Assuming that the above file is "input" and that our path is defined (boat is at /store/us/ron/moil.src/exe/boat ) we run boat < input > output The file "output" is below (self explanatory) boat> ~ BOAT TEST boat> ~ boat> ~debug boat> file rcon name=(VAL.WCON) unit=10 read boat> file rcrd name=(CMIN.CRD) unit=11 read boat> file boat name=(VAL.BOAT) unit=13 wovr boat> coor CHAR boat> action getcrd> * title for CHARMM coordinates getcrd> * **** BOAT COMPLETED **** The output file VAL.BOAT is below (self explanatory too) BOnd Angle Torsion (BOAT) list for structure # 1 For the bond, angle, improper torsion we provide the atomic indices, current value and deviation from equilibrium. Deviation not given for torsion **** BONDS 1 2 1.52150 -0.00050 2 4 1.33667 0.00167 2 3 1.22880 -0.00020 4 6 1.46210 0.01310 4 5 1.00675 -0.00325 6 7 1.55491 0.02891 6 10 1.53710 0.01510 7 8 1.53293 0.00693 7 9 1.52491 -0.00109 10 12 1.33179 -0.00321 10 11 1.22721 -0.00179 12 14 1.44775 -0.00125 12 13 1.02224 0.01224 **** ANGLES 1 2 3 118.78615 -1.61385 1 2 4 119.33635 -2.56365 3 2 4 121.74015 -1.15985 2 4 5 116.83170 -2.96830 2 4 6 128.11996 6.21996 5 4 6 115.03981 -3.36019 4 6 7 110.27675 0.57675 4 6 10 113.46651 3.36651 7 6 10 116.52774 5.42774 6 7 8 115.61151 4.11151 6 7 9 109.04338 -2.45662 8 7 9 110.49679 -1.00321 6 10 11 119.69924 -0.70076 6 10 12 117.46957 0.86957 11 10 12 122.58952 -0.31048 10 12 13 116.02447 -3.77553 10 12 14 124.12243 2.22243 13 12 14 119.82947 1.42947 **** TORSIONS 1 2 4 5 -1.26267 1 2 4 6 177.61388 3 2 4 5 174.42443 3 2 4 6 -6.69901 4 6 7 8 123.85246 4 6 10 11 129.39336 6 10 12 13 -3.83999 6 10 12 14 177.93716 11 10 12 13 -178.16854 11 10 12 14 3.60860 **** IMPROPER TORSIONS 2 1 3 4 -2.40613 -2.40613 4 6 2 5 -0.62396 -0.62396 6 4 10 7 28.82668 -6.43332 7 9 8 6 31.24737 -4.01263 10 6 11 12 3.22412 3.22412 12 14 10 13 0.97720 0.97720 CCRD.PRG ccrd is a moil program that conveniently(?!) allows you to change between format of coordinates allowed in moil. The are three type of data sets that moil recognize at present: (a) CHARMm formatted coordinate file (denoted by chr) (b) CHARMm and QUANTA dynamics format (unformatted single percision coordinated) - denoted by dyn and used primarily for display and analysis that does not invoke energy calculations (c) MOIL path format, double precision unformatted coordinates (denoted by pth) and are used in compact storage of coordinates when double precision is of importance, Typical input will look something like: ~ file conn name=(ALL.WCON) unit=10 read file rcrd name=(ALA11_20.PTH) binary unit=11 read file wcrd name=(20_5.DCD) unit=12 wovr fpth tchr lpst=5 action where we need to provide connectivity file (conn), file from which coordinates will be read (rcrd) and file to which the coordinates will be written. The line before action read in (...) we added interpretation of the keywords fpth (FROM PATH FORMAT) (WRITE) tchr (TO CHARMM) lpst=5 (LOOP START FROM 5) In other words pick structure number 5 in the path formatted file and write it in CHARMM format. Other combinations are possible, e.g. fdyn tpth lpst=2 lpend=51 in which structures 2 to 51 in the dynamics coordinate file are used to generate a new path file. Note that if dyn or pth files are used it is necessary to add the keyword "bina" in the file declaration. "Typical" output will look something like ccrd> ~ ccrd> file conn name=(ALL.WCON) unit=10 read ccrd> file rcrd name=(ALA11_20.PTH) binary unit=11 read ccrd> file wcrd name=(20_5.DCD) unit=12 wovr ccrd> fpth tchr lpst=5 ccrd> action lpst lpend 5 1 *** READING PATH FILE 5 ENERGY = -378.84814 CHAIN.PRG CHAIN is a program to find minimum energy paths by the method of simulated annealing. The functional to be optimized is similar to the SPW (Self Avoiding Walk) chain of Czerminski and Elber (Int. J. Quant. Chem. 24,167(1990)) except that only second nearest neighbor repulsion is taken into account. It is using a set of system copies interpolating between fixed reactant and product, the system copies are "intraconnected" via a set of equidistance constraints and repuslion between next nearest neighbours. c------------------------------------------------------------------------- Update: chmin is an alterantive program to chain that is less expensive that was also optimized considerably. The general structure is quite similar to chain but there are couple of important differences. It is explained in the tutorial for input sample see moil.tests/path c------------------------------------------------------------------------- Input files: (i) connectivity file , (ii) coordinate file for the complete path in path format (see below). All structures must be provided in the current version. In the near future a number of initiators will be added Output files: (i) standard output with information on the path progress (ii) "flashes" of the path coordinates during the simulated annealing progress in path format (hopefully in the near future CHARMm dynamics will be supported too. This is however less atractive for real calculation, since CHARMm dynamics file are only single precision) (iii) "flashes" of the path velocities, useful for restarting the job Path coordinate format : Ler r be the coordinate vector of the complete path. Let the total number of structures along the path be igrid and let a structure along the path be denoted by i (i=1,igrid). Furthermore the number of atoms is npt and an an atom is indexed by j, then r(1,...,npt) are the x coordinates of the first structure and r(npt+1,...,2*npt) - y coordinates and so on. A path is a sequence of igrid records of the follwoing type (all numbers are double precision) energy(i),r((i-1)*3*npt+1,..........,i*3*npt) where i=1,igrid and energy(i) is the energy of the i-th structure, kept for more convenient analysis In fortran, the following code will do the job k = -3*npt do i=1,igrid k = k + 3*npt write/read(upath)e0(i),(r(k+j),j=1,npt3) end do The following is a sample study of reaction path for conformational transition in valine dipeptide ~ CHAIN TEST FOR CONFORMATIONAL TRANSITION IN VALINE DIPEPTIDE ~ file conn name=(VAL.WCON) unit=10 read file rcrd name=(PATH.BIC) bina unit=11 read file wcrd name=(VALDYN.DCD) bina unit=12 wovr file wvel name=(VALDYN.VCD) bina unit=13 wovr #ste=10 #equ=1 #pri=1 #wcr=2 #lis=2 #sve=2 rand=-3451187 step=0.001 repl=10000. temp=100. lmbd=2. gama=200. rmax=9999. epsi=1. cdie v14f=8. el14=2. cpth action VAL.WCON is the connectivity file, PATH.BIC is the file with the path coordinates to start the calculation, VALDYN.DCD is a path formatted coordinate file and VALDYN.VCD, the same thing with velocities. The line below there is a mixture of energy parameters and chain parameters. #ste - number of dynamics steps #equ - number of equilibration steps (no velocity scaling) #pri - number of steps between printing out info #wcr - number of steps between writing coordinates #lis - number of steps between updates of the non-bonded list #sve - number of steps between velocity scaling rand - an initiator for velocity selection (Boltzmann distribution) repl - repulsion parameter between the structures (chain energy) temp - temperature of the simulation. A feature to be added is a continuous cooling. Currently it is necessary to restart at another temperature lmbd - lambda range parameter for second nearest neighbour repulsion gama - "spring" force constant for the equidistance constraint rmax - cutoff distance for nonbonded interactions epsi - dielectric constant v14f - 1-4 scaling factor for van der Waals interactions el14 - 1-4 scaling factor for electrostatic interactions cpth - coordinates to be written in PATH format (onlt way so far) action starts the action. Note that at the first step all the structures are oriented (minimum mass weighted rms) with respect to the middle structure. This is in order to effectively apply the rigid body constraints required in the CHAIN. See standard output file below. chain> ~ CHAIN TEST FOR CONFORMATIONAL TRANSITION IN VALINE DIPEPTIDE chain> ~ chain> file conn name=(VAL.WCON) unit=10 read chain> file rcrd name=(PATH.BIC) bina unit=11 read chain> file wcrd name=(VALDYN.DCD) bina unit=12 wovr chain> file wvel name=(VALDYN.VCD) bina unit=13 wovr chain> #ste=10 #equ=1 #pri=1 #wcr=2 #lis=2 chain> #sve=2 rand=-3451187 step=0.001 repl=10000. temp=100. chain> lmbd=2. gama=200. chain> rmax=9999. epsi=1. cdie v14f=8. e14f=2. cpth chain> action PARAMETERS FOR CHAIN SIMULATION: temperature: 100.0 number of integration steps: 10 print each 1 steps number of grid points of the chain : 21 data is on unit: 6 coordinates are on unit 12 initial coordinates are read from unit : 11 Coordinate are read in PATH style coordinates are written at step interval : 2 velocity scaling at step interval of: 2 step size: 0.100E-02 period for checking constraints 500 number of thermalization steps 1 select new velocities each 1000 steps initial seed for random vel. selection -3451187 Monomer - monomer parameters: i,i+1 force constant 200.00000 i,i+2 repulsion strength 10000.0000 i,i+2 repulsion range 2.00000 debug ? F rotation matrix 0.9131046947036799 -0.3945071968153701 -0.1029751823061233 0.3889723231938801 0.9185814987149735 -7.0061130505500668E-02 0.1222307174946419 2.3918651287898267E-02 0.9922134598065363 After overlap rms is = 1.077679313184675 rotation matrix 0.9281336179369756 -0.3610219971137616 -9.0725436649333357E-02 0.3586042785989994 0.9325334362964018 -4.2241704034229978E-02 9.9854687550024454E-02 6.6714158328969597E-03 0.9949796649103576 After overlap rms is = 1.029395180391329 rotation matrix 0.9414518729075730 -0.3282084956677283 -7.7120388810940105E-02 0.3280835326238476 0.9445351515068811 -1.4647292887457773E-02 7.7650284094009037E-02 -1.1512208276496745E-02 0.9969141901090142 After overlap rms is = 0.9473212813538049 rotation matrix 0.9533938794307759 -0.2939389845016896 -6.8117428416307552E-02 0.2950289126423913 0.9554675105520728 6.3068997513198652E-03 6.3230146045876706E-02 -2.6109570458746789E-02 0.9976573755360185 After overlap rms is = 0.8582695260021219 rotation matrix 0.9640816899876297 -0.2580999844540922 -6.2696834492610820E-02 0.2598962713614618 0.9653795055104770 2.2278655093070490E-02 5.4776118546368754E-02 -3.7773116963575638E-02 0.9977839287500326 After overlap rms is = 0.7648975293101804 rotation matrix 0.9734823007313952 -0.2211793636265584 -5.8411465214120953E-02 0.2232681618757305 0.9742326241242178 3.1970642543947304E-02 4.9835108661888838E-02 -4.4164275130365055E-02 0.9977805263417563 After overlap rms is = 0.6596051178409951 rotation matrix 0.9817587078525696 -0.1827294375036160 -5.2533724646796912E-02 0.1846102696558853 0.9822417564501140 3.3469093554745546E-02 4.5485029330918594E-02 -4.2556839114375501E-02 0.9980581283428944 After overlap rms is = 0.5453482660769873 rotation matrix 0.9889836341647168 -0.1413594746305001 -4.3918905798511709E-02 0.1426637775711453 0.9893716175156224 2.8122038042287731E-02 3.9476802346192497E-02 -3.4077872391186724E-02 0.9986392144762832 After overlap rms is = 0.4215237353632078 rotation matrix 0.9947993131182667 -9.6666983237275364E-02 -3.2090823785441258E-02 9.7296721499340821E-02 0.9950801845401103 1.8675500558253594E-02 3.0127638555049491E-02 -2.1700707072026860E-02 0.9993104646243257 After overlap rms is = 0.2879056201442189 rotation matrix 0.9986394659836653 -4.9226708965839194E-02 -1.7203142278409828E-02 4.9375215667922717E-02 0.9987456721207196 8.3168802971487002E-03 1.6772151251547851E-02 -9.1549737587569603E-03 0.9998174240319466 After overlap rms is = 0.1464982984040649 rotation matrix 1.000000000000000 -1.5764299587939234E-16 4.8572257327350599E-16 -2.2226144535952841E-16 1.000000000000000 -2.6020852139652106E-17 4.4408920985006262E-16 7.6327832942979512E-17 1.000000000000006 After overlap rms is = 0.0000000000000000E+00 rotation matrix 0.9985749006423962 5.0215098434739439E-02 1.8072401506679311E-02 -5.0141972631779420E-02 0.9987320648375199 -4.4771917514005079E-03 -1.8274309497846454E-02 3.5646254465731830E-03 0.9998266565049178 After overlap rms is = 0.1484263289405026 rotation matrix 0.9942182238547264 0.1013268110762974 3.5538721308762247E-02 -0.1012055164288685 0.9948519497184183 -5.2001524812332750E-03 -3.5882681052555113E-02 1.5733717203907892E-03 0.9993547706905281 After overlap rms is = 0.2952816549051642 rotation matrix 0.9866984354182601 0.1541295371091817 5.1674784311756712E-02 -0.1541532875188210 0.9880405679326983 -3.5496572187179928E-03 -5.1603890263204916E-02 -4.4633966595030265E-03 0.9986576573580991 After overlap rms is = 0.4358547145103406 rotation matrix 0.9754899602250702 0.2099319490232575 6.5938716088405305E-02 -0.2103087540215195 0.9776342156482107 -1.2523481243676936E-03 -6.4726852866542312E-02 -1.2645836200303687E-02 0.9978228887657316 After overlap rms is = 0.5619267869543590 rotation matrix 0.9599392773357913 0.2695939341427552 7.6391717492606342E-02 -0.2703470747082509 0.9627627998005902 -5.0051655476890389E-04 -7.3682040041861172E-02 -2.0171811856185905E-02 0.9970777577409439 After overlap rms is = 0.6558313416398356 rotation matrix 0.9402391552799410 0.3305108219943094 8.1932456469396453E-02 -0.3313541019955581 0.9435000632994414 -3.4770166330953742E-03 -7.8452469490646590E-02 -2.3879428355711149E-02 0.9966318191449750 After overlap rms is = 0.6984059235278118 rotation matrix 0.9153056842965772 0.3929597813985590 8.8306933462617643E-02 -0.3935730772873300 0.9192264354385485 -1.1090231085127404E-02 -8.5532082454222086E-02 -2.4604279996309875E-02 0.9960315719277727 After overlap rms is = 0.7266976123163730 rotation matrix 0.8823067159674221 0.4601228726960786 9.9104999776458855E-02 -0.4603774911787164 0.8874590950261144 -2.1654566987591398E-02 -9.7915394983499743E-02 -2.6519741275834618E-02 0.9948413334536761 After overlap rms is = 0.7663054209446161 rotation matrix 0.8392497322470157 0.5314220500442227 0.1151107799040020 -0.5316613974585185 0.8463826856385417 -3.1185058982544034E-02 -0.1140001990163503 -3.5027905705094634E-02 0.9928630320674348 After overlap rms is = 0.8112029601689894 rotation matrix 0.7882511892616290 0.6000976642279707 0.1361721558019867 -0.6013858664822667 0.7981414484571897 -3.6128490835245053E-02 -0.1303652646336490 -5.3413684040627976E-02 0.9900261997210925 After overlap rms is = 0.8445021186960484 current temperature is 106.072 chain energy at step 1 = -60.92350018991310 distances between monomers i & i+1 1 0.16006 2 0.16270 3 0.16109 4 0.15808 5 0.15635 6 0.15627 7 0.15551 8 0.15640 9 0.15930 10 0.16281 11 0.16547 12 0.16628 13 0.16456 14 0.15983 15 0.15171 16 0.14408 17 0.13905 18 0.13103 19 0.12110 20 0.11868 distances between monomers i & i+2 21 0.29855 22 0.31996 23 0.31622 24 0.30819 25 0.30809 26 0.30822 27 0.30759 28 0.31207 29 0.31934 30 0.32605 31 0.32973 32 0.32874 33 0.32146 34 0.30448 35 0.28014 36 0.27609 37 0.26786 38 0.25015 39 0.23409 chain energy at step 2 = -89.06535441169538 distances between monomers i & i+1 1 0.16629 2 0.16415 3 0.16169 4 0.15581 5 0.15450 6 0.15600 7 0.15386 8 0.16012 9 0.16155 10 0.16128 11 0.16661 12 0.16757 13 0.16399 14 0.15582 15 0.15549 16 0.14734 17 0.13743 18 0.12689 19 0.13665 20 0.13503 distances between monomers i & i+2 21 0.29876 22 0.32162 23 0.31300 24 0.30484 25 0.30540 26 0.30640 27 0.30919 28 0.31773 29 0.31942 30 0.32556 31 0.33187 32 0.32862 33 0.31681 34 0.30375 35 0.28502 36 0.27650 37 0.26110 38 0.25687 39 0.25583 step = 2 structure number - energy 1 -10.45415 2 -9.79856 3 -8.96012 4 -10.06318 5 -10.40488 6 -11.12723 7 -11.07365 8 -9.45259 9 -7.36740 10 -6.35559 11 -4.59451 12 -2.91222 13 -2.77387 14 -2.74422 15 -4.85994 16 -7.72668 17 -9.98231 18 -11.64531 19 -12.81461 20 -15.14641 21 -16.26571 *** scaling velocities at step 2 current temperature is 124.3667106044417 desired temperature is 100.0000000000000 scaling factor = 0.8967015616874476 chain energy at step 3 = -102.3447445399710 distances between monomers i & i+1 1 0.16740 2 0.16585 3 0.15996 4 0.15544 5 0.15418 6 0.15578 7 0.15348 8 0.16342 9 0.16369 10 0.16079 11 0.16766 12 0.16776 13 0.16333 14 0.15282 15 0.15811 16 0.15191 17 0.13600 18 0.12680 19 0.16799 20 0.12624 distances between monomers i & i+2 21 0.30244 22 0.32032 23 0.30958 24 0.30374 25 0.30409 26 0.30565 27 0.31176 28 0.32231 29 0.32042 30 0.32569 31 0.33233 32 0.32746 33 0.31194 34 0.30201 35 0.29060 36 0.27679 37 0.25760 38 0.27575 39 0.26866 chain energy at step 4 = -97.56480053161296 distances between monomers i & i+1 1 0.16843 2 0.16622 3 0.15647 4 0.15650 5 0.15488 6 0.15553 7 0.15447 8 0.16476 9 0.16453 10 0.16142 11 0.16791 12 0.16583 13 0.16190 14 0.15205 15 0.15684 16 0.15534 17 0.13812 18 0.13066 19 0.19182 20 0.12205 distances between monomers i & i+2 21 0.30549 22 0.31621 23 0.30674 24 0.30477 25 0.30410 26 0.30625 27 0.31421 28 0.32361 29 0.32164 30 0.32604 31 0.32984 32 0.32420 33 0.30910 34 0.29922 35 0.29272 36 0.28009 37 0.26167 38 0.29453 39 0.27645 step = 4 structure number - energy 1 -10.45415 2 -9.05355 3 -8.46168 4 -10.26709 5 -10.89047 6 -11.47964 7 -11.24178 8 -9.79647 9 -8.45058 10 -6.02135 11 -4.72257 12 -3.18838 13 -2.98129 14 -3.78811 15 -5.42577 16 -7.81969 17 -9.62130 18 -11.79652 19 -11.17097 20 -11.66508 21 -16.26571 chain energy at step 5 = -97.26216067862747 distances between monomers i & i+1 1 0.16842 2 0.16489 3 0.15308 4 0.15823 5 0.15602 6 0.15556 7 0.15658 8 0.16408 9 0.16421 10 0.16275 11 0.16714 12 0.16184 13 0.15996 14 0.15292 15 0.15256 16 0.15741 17 0.14494 18 0.13829 19 0.19501 20 0.12016 distances between monomers i & i+2 21 0.30644 22 0.31055 23 0.30467 24 0.30698 25 0.30503 26 0.30817 27 0.31578 28 0.32139 29 0.32239 30 0.32595 31 0.32413 32 0.31818 33 0.30794 34 0.29559 35 0.29143 36 0.28771 37 0.27420 38 0.30503 39 0.27626 chain energy at step 6 = -106.1717050497004 distances between monomers i & i+1 1 0.16714 2 0.16236 3 0.15174 4 0.15941 5 0.15723 6 0.15622 7 0.15887 8 0.16185 9 0.16281 10 0.16386 11 0.16523 12 0.15689 13 0.15787 14 0.15420 15 0.14911 16 0.15803 17 0.15468 18 0.14867 19 0.18036 20 0.12051 distances between monomers i & i+2 21 0.30434 22 0.30583 23 0.30382 24 0.30890 25 0.30649 26 0.31084 27 0.31557 28 0.31660 29 0.32167 30 0.32445 31 0.31650 32 0.31036 33 0.30684 34 0.29276 35 0.28978 36 0.29687 37 0.29284 38 0.30566 39 0.26996 step = 6 structure number - energy 1 -10.45415 2 -9.65276 3 -7.82695 4 -10.23818 5 -11.16747 6 -11.47818 7 -11.20173 8 -9.78400 9 -7.76895 10 -6.05645 11 -5.05099 12 -3.10875 13 -3.09628 14 -1.87249 15 -4.20945 16 -6.32040 17 -8.10151 18 -9.29663 19 -8.01742 20 -14.43783 21 -16.26571 chain energy at step 7 = -113.7698088467858 distances between monomers i & i+1 1 0.16449 2 0.15965 3 0.15292 4 0.15916 5 0.15832 6 0.15756 7 0.16018 8 0.15904 9 0.16041 10 0.16359 11 0.16232 12 0.15244 13 0.15559 14 0.15433 15 0.14970 16 0.15786 17 0.16364 18 0.15786 19 0.16209 20 0.12599 distances between monomers i & i+2 21 0.30001 22 0.30387 23 0.30416 24 0.30952 25 0.30813 26 0.31320 27 0.31345 28 0.31110 29 0.31850 30 0.32065 31 0.30892 32 0.30267 33 0.30405 34 0.29199 35 0.29143 36 0.30348 37 0.31080 38 0.30065 39 0.26379 chain energy at step 8 = -112.9418826856304 distances between monomers i & i+1 1 0.16014 2 0.15841 3 0.15590 4 0.15806 5 0.15915 6 0.15924 7 0.15973 8 0.15689 9 0.15785 10 0.16096 11 0.15902 12 0.14972 13 0.15367 14 0.15270 15 0.15426 16 0.15756 17 0.16954 18 0.16157 19 0.15124 20 0.13515 distances between monomers i & i+2 21 0.29552 22 0.30517 23 0.30525 24 0.30913 25 0.30975 26 0.31420 27 0.31019 28 0.30672 29 0.31286 30 0.31434 31 0.30285 32 0.29717 33 0.29976 34 0.29331 35 0.29677 36 0.30555 37 0.32139 38 0.29489 39 0.26371 step = 8 structure number - energy 1 -10.45415 2 -9.41156 3 -8.16995 4 -10.04024 5 -11.98986 6 -11.81939 7 -11.31188 8 -9.83658 9 -8.13657 10 -5.44244 11 -4.94568 12 -3.59278 13 -3.54205 14 -1.94643 15 -3.99327 16 -6.30790 17 -7.31663 18 -11.29136 19 -10.73587 20 -13.67689 21 -16.26571 chain energy at step 9 = -105.5633870812341 distances between monomers i & i+1 1 0.15431 2 0.16039 3 0.15864 4 0.15764 5 0.15967 6 0.16062 7 0.15769 8 0.15610 9 0.15623 10 0.15605 11 0.15611 12 0.14886 13 0.15308 14 0.15021 15 0.16005 16 0.15738 17 0.17224 18 0.15958 19 0.15184 20 0.14391 distances between monomers i & i+2 21 0.29273 22 0.30896 23 0.30657 24 0.30910 25 0.31118 26 0.31341 27 0.30709 28 0.30445 29 0.30605 30 0.30654 31 0.29879 32 0.29513 33 0.29623 34 0.29599 35 0.30296 36 0.30402 37 0.32308 38 0.29260 39 0.27184 chain energy at step 10 = -100.9488591626611 distances between monomers i & i+1 1 0.14839 2 0.16549 3 0.15905 4 0.15892 5 0.16003 6 0.16108 7 0.15513 8 0.15638 9 0.15562 10 0.15052 11 0.15403 12 0.14902 13 0.15412 14 0.14869 15 0.16365 16 0.15764 17 0.17225 18 0.15489 19 0.16184 20 0.15063 distances between monomers i & i+2 21 0.29226 22 0.31354 23 0.30796 24 0.31071 25 0.31221 26 0.31121 27 0.30506 28 0.30406 29 0.29993 30 0.29922 31 0.29649 32 0.29650 33 0.29592 34 0.29918 35 0.30682 36 0.30134 37 0.31878 38 0.29564 39 0.28527 step = 10 structure number - energy 1 -10.45415 2 -10.86950 3 -8.34173 4 -9.47768 5 -11.30336 6 -10.78043 7 -10.64849 8 -9.18179 9 -6.76493 10 -5.48701 11 -4.81730 12 -3.04506 13 -2.72662 14 0.28894 15 -1.85944 16 -5.27395 17 -6.47830 18 -9.20206 19 -4.67340 20 -15.18820 21 -16.26571 CON.PRG con is a program to form molecules from a data base. I.e. it assigns particle properties bonds, etc according to the audience request (if everything works as planned). It uses three data sets: (i) the property file, (ii) the monomer file and (iii) the polymerization file (see also prop.file, mono.file, and poly.file in the present directory) to kindly provide the user with a connectivity file (see conn.file). Optionally, it can also employ two more data sets: addbond file (unit ubond) and edcon (uedi). addbond can be used to add bond(s) to an existing structure. Useful for creating S-S bonds and linkages between groups like heme and protein. The syntax of addbond is bond atm1=[int] atm2=[int] where [int] is an integer that is the particle index in the final structure. The first number (atm1) MUST be smaller than the second number (atm2). EXAMPLE: bond atm1=10 atm2=100 OR bond chem [mono-name] [mono-index] [unique-particle] - [mono-name] [mono-index] [unique-particle] chem - meaning bond according to chemical definition of particles mono-name - monomer name used to identify the particle mono-index - number of the monomer to be bonded unique particle - unique name of particle in the given monomer The above triplet is used to find the particle absolute number the last are repeated to identify the second particle to be bonded to. EXAMPLE: bond chem HIS 95 NE2 HEM1 157 FE (make a bond between Ne2 of His 95 to the iron atom of the Hem group) The second option is edcon which suppose to make it possible to modify the structure after creation. Currently it only allows to remove angles (a response to an ergent need). Future expansion expected. Syntax remo angl atm1=[int] atm2=[int] atm3=[int] Below breif description of what is in the different files is given. For more details consults the above mentioned *.file. (i) property file. Must be read first. Includes definition of particle types as mass, charge, van der Waals radius etc. and covalent structure parameters as bond length or improper torsion force constants. (ii) monomer file. includes definition of monomers to be linked together to one big happy molecule. For each monomer unique particles are provided and also the set of bond inside that monomer and to its neighbors. (iii) polymerization file. where the sequence of the actual monomers to be used is provided. with those three what else one can ask? the output of the connectivity in the molecule and the non-bonded interactions are written on the connectivity file. An input file to con that basically direct con of where to read and write those files is given below ~ Input to generate valine dipeptide ~ (If you uncomment the line below you will fill the disk with ~ debugging information) ~debug file prop name=(ALL.PROP) unit=10 read file mono name=(VAL.MONO) unit=11 read file poly name=(VAL.POLY) unit=15 read file wcon name=(VAL.WCON) unit=13 wovr ~ optional files are commented below ~ file ubon name=(whatever) unit=16 read ~ file uedi name=(never) unit=17 read action *EOD the "file" commands give direction where to find the differet files. For example "file prop" refers to the property file. The program will find it at ALL.PROP, will open it for read only and will attach it internally to unit 10. You may use any units you like, excluding unit 25 that is employed in the interpretation of the command line. "wovr" means that if the file already exists write over it. Other option is "write". This option bounces if you try to read on an existing file. "wcon" as you should have guess by now is the connectivity file. After all the file assingments end the command "action", initiate somore serious activity hopefully resulting in an appropriate connecntivity file. This code is doing a lot of i/o in line interpretation, so some of the stuff may look quite slow. However, the program should terminate after a few seconds for reasonably sized molecule. Set your path variable to include (on bully) /store/us/ron/moil.src/exe To run the above input type (assuming it is in file "input") con < input > output The output file that also lists what was in the different input data base is given below CONN> ~ Input to generate valine dipeptide CONN> ~ (If you uncomment the line below you will fill the disk with CONN> ~ debugging information) CONN> ~debug CONN> file prop name=(ALL.PROP) unit=10 read CONN> file mono name=(VAL.MONO) unit=11 read CONN> file poly name=(VAL.POLY) unit=15 read CONN> file wcon name=(VAL.WCON) unit=13 wovr CONN> ~ optional files are commented below CONN> ~ file ubon name=(whatever) unit=16 read CONN> ~ file uedi name=(never) unit=17 read CONN> action propert> PRTC propert> ~ propert> ~ The info below is taken from OPLS paper propert> ~ Jorgensen and Tirado-Rives, JACS 110,1657(1988) propert> ~ propert> ~ structure of data propert> ~ name mass charge epsilon sigma propert> ~ propert> ~ CANX NX & HX are N-terminal particles propert> ~ propert> ~ HX propert> ~ / propert> ~ CANX-NX-HX propert> ~ | \ propert> ~ CB HX propert> ~ propert> PNAM=(NX) PMAS=14. PCHG=-0.3 PEPS=0.170 PSGM=3.250 propert> PNAM=(HX) PMAS=1. PCHG=0.33 PEPS=0.0 PSGM=0.0 propert> PNAM=(CANX) PMAS=13. PCHG=0.310 PEPS=0.118 PSGM=3.800 propert> ~ propert> ~ CACX CX & OX are C-terminal particles propert> ~ propert> ~ OX propert> ~ / propert> ~ CACX-CX propert> ~ | \ propert> ~ CB OX propert> ~ propert> PNAM=(CX) PMAS=12. PCHG=0.70 PEPS=0.105 PSGM=3.750 propert> PNAM=(OX) PMAS=16. PCHG=-0.80 PEPS=0.210 PSGM=2.960 propert> PNAM=(CACX) PMAS=13. PCHG=0.10 PEPS=3.80 PSGM=0.080 propert> ~ propert> ~ Normal backbone includes: NH HN CO OC & CAH propert> ~ propert> ~ HN OC propert> ~ | | propert> ~ NH-CAH-CO propert> ~ | propert> ~ CB propert> ~ propert> PNAM=(NH) PMAS=14. PCHG=-0.570 PEPS=0.105 PSGM=3.250 propert> PNAM=(HN) PMAS=1. PCHG=0.370 PEPS=0.00 PSGM=0.00 propert> PNAM=(CO) PMAS=12. PCHG=0.50 PEPS=0.105 PSGM=3.750 propert> PNAM=(OC) PMAS=16. PCHG=-0.5 PEPS=0.21 PSGM=2.960 propert> PNAM=(CAH) PMAS=13. PCHG=0.20 PEPS=0.080 PSGM=3.800 propert> ~ propert> ~ CH3 is as suggested propert> ~ propert> PNAM=(CH3) PMAS=15. PCHG=0.00 PEPS=0.160 PSGM=3.910 propert> ~ propert> ~ CH3G - dummy atom for valine monomer propert> ~ propert> PNAM=(CH3G) PMAS=15. PCHG=0.00 PEPS=0.160 PSGM=3.910 propert> ~ CBH - B carbon for branched side chain (e.g. valine) propert> ~ propert> PNAM=(CBH) PMAS=13. PCHG=0.00 PEPS=0.080 PSGM=3.850 propert> ~ propert> ~ CH3T used for neutral terminus from the nitrogen side propert> ~ propert> PNAM=(CH3T) PMAS=15. PCHG=0.200 PEPS=0.170 PSGM=3.800 propert> DONE *** DONE WITH PARTICLE PROPERTIES *** propert> ~ propert> ~ The covalent part of this potential is taken from AMBER paper propert> ~ Weiner et al. JACS 106,765(1984) propert> ~ The name of the particles is however different. propert> ~ propert> BOND propert> NX HX 434.0 1.01 propert> CANX NX 337.0 1.449 propert> CANX CH3 260.0 1.526 propert> CANX CO 317.0 1.522 propert> CACX NH 337.0 1.449 propert> CACX CH3 260.0 1.526 propert> CACX CX 317.0 1.522 propert> CX OX 656.0 1.250 propert> CAH CH3 260.0 1.526 propert> CAH CO 317.0 1.522 propert> CO OC 570.0 1.229 propert> CAH CBH 260.0 1.526 propert> CAH NH 337.0 1.449 propert> NH HN 434.0 1.010 propert> CO CH3 317.0 1.522 propert> CO NH 490.0 1.335 propert> CBH CH3 260.0 1.526 propert> CBH CH3G 260.0 1.526 propert> CH3T NH 337.0 1.449 propert> DONE *** DONE WITH BOND PROPERTIES *** propert> ~ propert> ANGLE propert> ~ propert> ~ According to AMBER, there seem to be no difference between propert> ~ angles of the type CH1-X-X, CH2-X-X and CH3-X-X . Somewhat propert> Strange... propert> ~ propert> HX NX HX 35.0 109.5 propert> HX NX CANX 35.0 109.5 propert> NX CANX CH3 35.0 109.5 propert> NX CANX CO 63.0 110.1 propert> CH3 CAH CO 63.0 111.1 propert> CH3 CANX CO 63.0 111.1 propert> CANX CO OC 80.0 120.4 propert> CANX CO NH 70.0 116.6 propert> CO NH CACX 50.0 121.9 propert> CACX NH HN 38.0 118.4 propert> NH CAH CH3 80.0 109.5 propert> NH CACX CH3 80.0 109.5 propert> NH CACX CX 63.0 110.1 propert> CH3 CACX CX 63.0 111.1 propert> CACX CX OX 65.0 117.0 propert> OX CX OX 80.0 126.0 propert> CH3 CO OC 80.0 120.4 propert> OC CO NH 80.0 122.9 propert> CH3 CO NH 50.0 121.9 propert> CO NH HN 35.0 119.8 propert> CO NH CAH 50.0 121.9 propert> CAH NH HN 38.0 118.4 propert> NH CAH CBH 80.0 109.7 propert> NH CAH CO 63.0 110.1 propert> CAH CBH CH3 63.0 111.5 propert> CAH CBH CH3G 63.0 111.5 propert> CH3 CBH CH3G 63.0 111.5 propert> CBH CAH CO 63.0 111.1 propert> CAH CO OC 80.0 120.4 propert> CO NH CH3T 50.0 121.9 propert> CH3T NH HN 38.0 118.4 propert> CAH CO NH 70.0 116.6 propert> DONE *** DONE WITH ANGLE PROPERTIES *** propert> ~ propert> TORSION propert> ~ propert> ~ Note that generic torsions as X-C-C-X must come at the end of propert> ~ torsion definition, since the program pick up the first compatible propert> ~ torsion. So generic definition must come last. propert> ~ propert> ~ The torsion parameters are derived from AMBER using different propert> ~ notation. I.e. we write E = SUM V(n)/2( 1 + cos(n*phi+gamma) propert> ~ as k(n)=V(n)/2 and phase=cos(gamma) propert> ~ Our data structure is therefore propert> ~ name name name name k(1) k(2) k(3) n cos(gamma) propert> ~ should be added cos(g2) cos(g3) propert> ~ propert> CAH CO NH CAH 0.0000 2.5000 0.0000 2 -1.0 propert> CO NH CAH CH3 0.0 0.0 0.0 3 0.0 propert> HN NH CAH CH3 0.0 0.0 0.0 3 0.0 propert> CANX CO NH CAH 0.0000 2.5000 0.0000 2 -1.0 propert> CH3 CO NH HN 0.0000 2.5000 0.0000 2 -1.0 propert> OC CO NH HN 0.0000 2.5000 0.0000 2 -1.0 propert> CH3 CO NH CAH 0.0000 2.5000 0.0000 2 -1.0 propert> OC CO NH CAH 0.0000 2.5000 0.0000 2 -1.0 propert> OC CO NH CH3T 0.0000 2.5000 0.0000 2 -1.0 propert> CAH CO NH CH3T 0.0000 2.5000 0.0000 2 -1.0 propert> CAH CO NH HN 0.0000 2.5000 0.0000 2 -1.0 propert> OC CO CAH NH 0.0 0.0 0.1 3 -1.0 propert> CO NH CAH CBH 0.0 0.0 0.0 3 0.0 propert> HN NH CAH CBH 0.0 0.0 0.0 3 0.0 propert> CO NH CAH CO 0.0 0.0 0.0 3 0.0 propert> HN NH CAH CO 0.0 0.0 0.0 3 0.0 propert> NH CAH CBH CH3G 0.0 0.0 2.0 3 -1.0 propert> CAH CO NH CACX 0.0000 2.5000 0.0000 2 -1.0 propert> CO NH CACX CH3 0.0 0.0 0.0 3 0.0 propert> CO NH CACX CX 0.0 0.0 0.0 3 0.0 propert> HN NH CACX CX 0.0 0.0 0.0 3 0.0 propert> CH3 CACX CX OX 0.0 0.0 1.0 3 1.0 propert> NH CACX CX OX 0.0 0.0 1.0 3 1.0 propert> CH3 CANX CO NH 0.0 0.0 0.0 3 0.0 propert> CH3 CANX CO OC 0.0 0.0 0.0 3 0.0 propert> OC CO NH CACX 0.0000 2.5000 0.0000 2 -1.0 propert> CANX CO NH HN 0.0000 2.5000 0.0000 2 -1.0 propert> CANX CO NH CACX 0.0000 2.5000 0.0000 2 -1.0 propert> HN NH CACX CH3 0.0 0.0 0.0 3 0.0 propert> HX NX CANX CH3 0.0 0.0 0.0 3 0.0 propert> HX NX CANX CO 0.0 0.0 0.0 3 0.0 propert> NX CANX CO OC 0.0 0.0 0.1 3 -1.0 propert> NX CANX CO NH 0.0 0.0 0.0 3 0.0 propert> X CAH CBH X 0.0 0.0 0.0 3 0.0 propert> X CAH CO X 0.0 0.0 0.0 3 0.0 propert> X CANX CX X 0.0 0.0 0.0 3 0.0 propert> DONE *** DONE WITH TORSION PROPERTIES *** propert> ~ propert> IMPROPER propert> ~ the improper are set to very high value since propert> ~ the current setup is such v(phi) = k (cos(phi) - cos(phieq))^2 propert> ~ which is computationally more convenient. propert> ~ yielding in the neighborhood of phieq propert> ~ k (phi+phieq)^2(phi-phieq)^2 . For phieq=0 and close to phieq propert> ~ this term is of fourth order propert> ~ propert> NH CAH CO HN 4500.0 0.0 propert> NH CH3T CO HN 4500.0 0.0 propert> CO CAH OC NH 10000.0 0.0 propert> CO CH3 OC NH 10000.0 0.0 propert> CO CANX OC NH 10000.0 0.0 propert> CX CACX OX OX 10000.0 0.0 propert> CX OX CACX OX 10000.0 0.0 propert> CANX NX CO CH3 55.0 35.26 propert> CACX NH CX CH3 55.0 35.26 propert> CAH NH CO CBH 55.0 35.26 propert> CAH NH CO CH3 55.0 35.26 propert> CBH CH3 CH3G CAH 55.0 35.26 propert> DONE *** DONE WITH IMPROPER TORSIONS *** *** Number of parameters read: Unique particles, bonds, angles, torsions, improper torsions 15 19 32 36 12 monomers> MONO LIST monomers> ~ monomers> MONO=(VALD) #prt=14 chrg=0 monomers> ~ monomers> ~ H1 O2 monomers> ~ | | monomers> ~ ME1-C1-N1-CA-C2-N2-ME2 monomers> ~ | | | monomers> ~ O1 CB H2 monomers> ~ / \ monomers> ~ CG1 CG2 monomers> ~ monomers> UNIQ=(ME1) PRTC=(CH3) monomers> UNIQ=(C1) PRTC=(CO) monomers> UNIQ=(O1) PRTC=(OC) monomers> UNIQ=(N1) PRTC=(NH) monomers> UNIQ=(H1) PRTC=(HN) monomers> UNIQ=(CA) PRTC=(CAH) monomers> UNIQ=(CB) PRTC=(CBH) monomers> UNIQ=(CG1) PRTC=(CH3G) monomers> UNIQ=(CG2) PRTC=(CH3) monomers> UNIQ=(C2) PRTC=(CO) monomers> UNIQ=(O2) PRTC=(OC) monomers> UNIQ=(N2) PRTC=(NH) monomers> UNIQ=(H2) PRTC=(HN) monomers> UNIQ=(ME2) PRTC=(CH3T) monomers> DONE monomers> BOND monomers> ME1-C1 C1-O1 C1-N1 N1-H1 N1-CA CA-CB CB-CG1 CB-CG2 monomers> CA-C2 C2-O2 C2-N2 N2-H2 N2-ME2 monomers> DONE monomers> *EOD LEGO> ~ LEGO> MOLC=(VALD) #mon=1 LEGO> VALD *** MOLECULE VALD ASSEMBLED *** LEGO> *EOD Note that the different subsection of the programs are echoing their names as things pass alone. Once the message "MOLECULE XXX ASSEMBLED" appears you are porbably ok. The careful user will examine the resulting connectivity file - VAL.WCON - to see if it makes sense. DYNA.PRG The dynamics program is using similar connectivity and coordinate files (coordinates are in CHARMm format) to run normal Verlet type MD. Output files are standard output that contains info on the simulation progress, velocity and coordinate files (CHARMm dynamics format). Below is a sample input ~ DYNAMICS TEST ~ file conn name=(VAL.WCON) unit=10 read file rcrd name=(A.CRD) unit=11 read file wcrd name=(VALDYN.DCD) bina unit=12 wovr file wvel name=(VALDYN.VCD) bina unit=13 wovr #ste=100 #equ=10 info=1 #crd=10 #vel=10 #lis=10 #scl=10 #tmp=1 rand=-3451187 step=0.001 rmax=9999. epsi=1. cdie v14f=8. e14f=2. shkb action The "bina" keyword in the file lines means that the file should be opened unformatted. Parameters that come after the file lines include normal energy parameters (e.g. rmax epsi) and the file energy.prog should be consulted for those. In addition simulation parameters include #ste - number of steps #equ - number of equilibration steps info - number of steps between writing information on standard output #crd - number of steps between writing coordinate sets #vel - number of steps between writing velocities sets #lis - number of steps between regeneration of the nonbonded list #scl - number of steps between velcoity scaling for requested temperature(s) #tmp - number of temperatures in the system (useful for LES simulation) if larger than one more input is required to define different domains with different temperatures (see below) tmpi - initial temperature tmpf - final temperature For a number of tmpi and tmpf for #tmp>1 rand - an initiator for (Boltzmann) velocity generator. Number must be negative step - time step in ps shkb - shake all bonds (alternatively one may try shkl for shaking bonds with light particles only m<1.1, shkb is highly recommended for dynamics. Also supported is particle selection. The selection will be used for multiple temperatures and "freezing" part of the molecule. nfrz - select the particles that WILL NOT be freezed. In the same line as nfrz a "pick" command must follow (see pick.sub and an example below for temperature selection) The line with "action" is self explanatory For LES run with multiple copies (as defined in the wcon generation) it is useful to define multiple temperature done via velocity scaling. Below is an input sample ~ ~ DYNAMICS TEST ~ file conn name=(wcon.mb.multco) unit=10 read file rcrd name=(MBMULTCO.MIN.CRD) unit=11 read file wcrd name=(MBMULTT.DCD) bina unit=12 wovr file wvel name=(MBMULTT.VCD) bina unit=13 wovr #ste=5000 #equ=500 info=50 #crd=20 #vel=50 #lis=10 #scl=10 #tmp=2 rand=-3451187 step=0.002 rmax=10. epsi=1. cdie v14f=8. e14f=2. shkb mult pick grou 1 chem prtc * | - grou 2 chem mono CO done - tmpi 300. 30. tmpf 300. 30. action where the new fautures start with "mult", note that in this case the program expect two temperatures (#tmp=2) and that the "mult" is a single command continuation using "-" is used. The structure is mult (state that multiple temperature command is coming) pick grou 1 (pick group 1 that belongs to temperature 1 end first pick by an or (|) sign grou 2 (pick particles that belong to group two for temperature assignment. See pick.sub for specific pick selection tmpi (initial temperatures of all the group listed in order) tmpf (final temperatures of all the groups listed in order) *** Note that in LES the actual equilibrium temperature of the *** enhanced part is N times higher N number of copies *** i.e., for 10 copies the enhanced temperatures should be set *** to 10 times smaller than the bath Below is a sample dyna output obtained from the benchmark run dynamics> ~ DYNAMICS TEST dynamics> ~ dynamics> file conn name=(wcon.new3.les5) unit=10 read dynamics> file rcrd name=(LEG5CO.CRD) unit=11 read getcrd> * title for CHARMM coordinates getcrd> * dynamics> file wcrd name=(LEG5CO.DCD) bina unit=12 wovr dynamics> file wvel name=(LEG5CO.VCD) bina unit=13 wovr dynamics> #ste=500 #equ=100 info=100 #crd=10000 #vel=10000 #lis=25 dynamics> #scl=25 #tem=1 rand=-3451187 step=0.003 tmpi=10 tmpf=300 dynamics> rmax=10. epsi=1. cdie v14f=8. e14f=2. shkb All bonds will be shaked Number of shake constraints = 1570 dynamics> action Number of degf (exc. shake) 2891 Parameter for energy calculation Constant dielectric will be used. Cutoff= 10.00000 Dielectric constant= 1.00000 1-4 scaling factors van der Waals= 8.00000 Electrostatic= 2.00000 Number of particle pairs 194174 Number of particle pairs 194959 Number of particle pairs 195606 Number of particle pairs 195072 At dynamics step 100 Current temperature(s) are 291.51 current energy (kinetic+potential) is -6555.064 ENERGIES: E total = -7392.675 E bond = .000 E angl = 517.011 E tors = 401.779 E impr = 118.263 E vdw = -562.139 E elec = -7867.591 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = 9.196 Number of particle pairs 193947 rotation matrix .9999851890726013 -4.297308194279601E-03 -3.339876901658505E-03 4.310177090564005E-03 .9999832785834191 3.855506063320246E-03 3.323252756386528E-03 -3.869844420607005E-03 .9999869900630097 After overlap rms is = 1.02984093077199 Scaling velocities at current temperatures 277.0994 Number of particle pairs 191872 rotation matrix .9999999938076689 -5.954827930246243E-07 -1.112847976718574E-04 6.127585633186182E-07 .9999999879501933 1.552392929161650E-04 1.112847038883401E-04 -1.552393601455537E-04 .9999999817582277 After overlap rms is = 1.19055593517002 Scaling velocities at current temperatures 284.0347 Number of particle pairs 188241 rotation matrix .9999999923305885 3.301705106196827E-05 -1.193678987038706E-04 -3.301478896786049E-05 .9999999992754127 1.895252617672982E-05 1.193685243740735E-04 -1.894858512435648E-05 .9999999926960531 After overlap rms is = 1.34913710999694 Scaling velocities at current temperatures 281.4057 Number of particle pairs 184520 At dynamics step 200 Current temperature(s) are 283.74 current energy (kinetic+potential) is -7071.470 ENERGIES: E total = -7886.748 E bond = .000 E angl = 543.578 E tors = 393.048 E impr = 124.240 E vdw = -567.911 E elec = -8379.703 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = 9.322 rotation matrix .9999999986285075 5.592775995100164E-06 -5.207403865474590E-05 -5.597945620458870E-06 .999999995056577 -9.927491263389542E-05 5.207348317484901E-05 9.927520400593303E-05 .9999999937163922 After overlap rms is = 1.49825295170669 Scaling velocities at current temperatures 283.7363 Number of particle pairs 180736 rotation matrix .9999999976850062 5.784195817772470E-05 3.583705259539791E-05 -5.783377460459007E-05 .9999999722639771 -2.283140361908353E-04 -3.585025773248317E-05 2.283119630711039E-04 .9999999732942037 After overlap rms is = 1.64159339243577 Scaling velocities at current temperatures 295.8566 Number of particle pairs 177401 rotation matrix .9999999995356015 2.793826804650678E-05 1.217576617833487E-05 -2.793544295454737E-05 .9999999727061158 -2.319641746007516E-04 -1.218224652360221E-05 2.319638343560859E-04 .9999999730221856 After overlap rms is = 1.77742081657301 Scaling velocities at current temperatures 267.3479 Number of particle pairs 174651 rotation matrix .9999999991042222 -2.796954405171977E-05 3.176883476161806E-05 2.798229401394835E-05 .9999999190450938 -4.014060223496707E-04 -3.175760504627468E-05 4.014069109545825E-04 .9999999189319692 After overlap rms is = 1.9038470956868 Scaling velocities at current temperatures 286.8885 Number of particle pairs 172571 At dynamics step 300 Current temperature(s) are 287.20 current energy (kinetic+potential) is -6596.078 ENERGIES: E total = -7421.316 E bond = .000 E angl = 571.700 E tors = 364.204 E impr = 120.739 E vdw = -465.740 E elec = -8012.220 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = 9.932 rotation matrix .9999999969829339 -3.234610882297010E-05 7.062478944903742E-05 3.237654922598302E-05 .999999906570699 -4.310572496377118E-04 -7.061083982592888E-05 4.310595349245049E-04 .9999999046008894 After overlap rms is = 2.01460169909338 Scaling velocities at current temperatures 287.2027 Number of particle pairs 171525 rotation matrix .9999999874039417 2.854105135304809E-05 1.561330312686466E-04 -2.848533799330518E-05 .9999999359328875 -3.568232133472349E-04 -1.561432053754138E-04 3.568187613507411E-04 .9999999241498323 After overlap rms is = 2.10585797679896 Scaling velocities at current temperatures 296.4653 Number of particle pairs 170510 rotation matrix .99999997932775 9.776958248996570E-05 1.782851885431171E-04 -9.771728402142247E-05 .9999999522053716 -2.933267532286287E-04 -1.783138584565509E-04 2.933093256197172E-04 .9999999410869007 After overlap rms is = 2.18103218658154 Scaling velocities at current temperatures 290.9355 Number of particle pairs 169690 rotation matrix .9999999567654868 1.125865227369123E-04 2.716492206370855E-04 -1.125535440459806E-04 .9999999862951727 -1.214139816524961E-04 -2.716628864923472E-04 1.213834013205628E-04 .9999999557326713 After overlap rms is = 2.24894443487863 Scaling velocities at current temperatures 292.7568 Number of particle pairs 169293 At dynamics step 400 Current temperature(s) are 293.01 current energy (kinetic+potential) is -6837.544 ENERGIES: E total = -7679.465 E bond = .000 E angl = 598.800 E tors = 419.164 E impr = 111.797 E vdw = -449.226 E elec = -8359.998 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = 10.038 rotation matrix .9999999468380948 1.688104565545898E-04 2.789746209579200E-04 -1.687638810240681E-04 .9999999718202726 -1.669676816704046E-04 -2.790027989867417E-04 1.669205919536675E-04 .9999999471474757 After overlap rms is = 2.31115352218454 Scaling velocities at current temperatures 293.0085 Number of particle pairs 169057 rotation matrix .9999999312811581 2.663847506553626E-04 2.578310387584097E-04 -2.663531276514306E-04 .9999999570032214 -1.226766702328574E-04 -2.578637068669943E-04 1.226079876981422E-04 .9999999592367937 After overlap rms is = 2.36394648837109 Scaling velocities at current temperatures 291.0874 Number of particle pairs 168407 rotation matrix .9999999660565705 2.327283407500325E-04 1.171510906354013E-04 -2.327221690037319E-04 .999999971531929 -5.269281825981409E-05 -1.171633504120431E-04 5.266555281596696E-05 .9999999917495451 After overlap rms is = 2.40857301362181 Scaling velocities at current temperatures 296.4393 Number of particle pairs 168026 rotation matrix .999999963589042 2.318163825954810E-04 1.381415175569691E-04 -2.318139400088147E-04 .9999999729745471 -1.769752248659273E-05 -1.381456163994289E-04 1.766549871312639E-05 .9999999903018599 After overlap rms is = 2.43912280500919 Scaling velocities at current temperatures 289.7260 Number of particle pairs 167564 At dynamics step 500 Current temperature(s) are 292.41 current energy (kinetic+potential) is -6682.757 ENERGIES: E total = -7522.955 E bond = .000 E angl = 583.521 E tors = 385.462 E impr = 118.360 E vdw = -447.453 E elec = -8162.844 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = 10.865 rotation matrix .9999999433104272 3.112712686524051E-04 1.284108235925715E-04 -3.112973517396533E-04 .9999999309115642 2.031522268078900E-04 -1.283475792693104E-04 -2.031921892398091E-04 .9999999711199171 After overlap rms is = 2.46749855422739 Scaling velocities at current temperatures 292.4090 real 12:13.33 user 12:12.68 sys 0.20 MINI.PRG The minimization program is similar in style to the energy program (recommended to read the file energy.prog first). It is used as input the standard input and two other files: (i) the connectivity file and (ii) coordinate file (currently only CHARMm format is supported) output files include the standard output, coordinates of minimized structure, and minimization history. "Typical" input is below ~ MINIMIZATION TEST ~ file rcon name=(ALL.WCON) unit=10 read file rcrd name=(DIALA.CRD) unit=11 read file wcrd name=(ALAMIN.CRD) unit=12 wovr file wmin name=(DIALA.MIN) unit=13 wovr rmax=9999. epsi=50. cdie v14f=8. e14f=2. tolf=0.01 mistep=200 list=50 action The file ALL.WCON is the connectivity file, DIALA.CRD is the initial coordinate set. The final coordinate set is ALAMIN.CRD and the minimization history is in DIALA.MIN. The line after gives energy parameters: rmax - cutoff distance, epsi - dielectric constant, cdie - use constant dielectric (versus distance dependent dielectric), v14f - scaling factor for 1-4 van der Waals interactions, e14f - scaling factor for 1-4 electrostatic interactions. More parameters are given in the next line (these parameters can be mixed with the energy parameters) , tolf - the requested gradient rms, mistep - number of minimization steps, list - number of steps between updates on the nonbonded list. We can also turn off diferent energy terms (see energy.prog) Running the minimization routine from input file "input" to output file "output" with a path defined properly to /store/us/ron/moil/src/exe is mini < input > output The output file using the above input is below mini> ~ MINIMIZATION TEST mini> ~ mini> ~debug mini> file rcon name=(ALL.WCON) unit=10 read mini> file rcrd name=(DIALA.CRD) unit=11 read mini> file wcrd name=(ALAMIN.CRD) unit=12 wovr mini> file wmin name=(DIALA.MIN) unit=13 wovr mini> rmax=9999. epsi=50. cdie v14f=8. e14f=2. mini> tolf=0.01 mistep=200 list=50 mini> action getcrd> * COORDINATES FOR ALANINE DIPEPTIDE getcrd> * DATE: 23/12/1991 11:42:43 CREATED BY USER: ron getcrd> * *** Minimization completed The DIALA.MIN file is below At step 0 Parameter for energy calculation Constant dielectric will be used. Cutoff= 9999.00000 Dielectric constant= 50.00000 1-4 scaling factors van der Waals= 8.00000 Electrostatic= 2.00000 After 50 energy calls Current gradient norm 2.08663 requested 0.01000 ENERGIES: E total = 2.09171 E bond = 0.03909 E angl = 0.43055 E tors = 3.26284 E impr = 0.16778 E vdw = 0.27732 E elec = -2.08587 Norm Force = 2.08663 After 100 energy calls Current gradient norm 0.81426 requested 0.01000 ENERGIES: E total = 1.19510 E bond = 0.01327 E angl = 0.03398 E tors = 3.03439 E impr = 0.25563 E vdw = -0.12398 E elec = -2.01818 Norm Force = 0.81426 After 150 energy calls Current gradient norm 0.64623 requested 0.01000 ENERGIES: E total = 1.02379 E bond = 0.01671 E angl = 0.05973 E tors = 3.10110 E impr = 0.22976 E vdw = -0.38536 E elec = -1.99815 Norm Force = 0.64623 After 200 energy calls Current gradient norm 0.36390 requested 0.01000 ENERGIES: E total = 0.95400 E bond = 0.00604 E angl = 0.02138 E tors = 3.08319 E impr = 0.24736 E vdw = -0.41422 E elec = -1.98975 Norm Force = 0.36390 RMS.PRG Here is how to use the rms program to calculate rms and fluctuation Needed file are: connectivity , coordintes of the reference structure and coordinates obtained during dynamics simulation . below is a sample input ~ ~ RMS TEST ~ file conn name=(wcon.mbco) unit=10 read file rcrd name=(MBCO.MIN.CRD) unit=11 read file rdyc name=(MBCO.DCD) bina unit=12 read file wrms name=(rms.data) unit=14 wovr file wave name=(rms_ave.data) unit=15 wovr file wflu name=(flu.data) unit=16 wovr #crd=100 #ste=1000 step=0.002 action #ste - number of steps in the dynamics simulation #crd - number of steps between coordinate sets step - time step in ps rms.data - result rms between dynamics structure and reference structure rms_ave.data - result rms between dynamics structrue and average structure flu.data - fluctuation from average struction ( (r-)^2 ) If you only want one of them , what you need to do is remove other two lines in the input. For example,if you want only rms.data, the input is ~ ~ RMS TEST ~ file conn name=(wcon.mbco) unit=10 read file rcrd name=(MBCO.MIN.CRD) unit=11 read file rdyc name=(MBCO.DCD) bina unit=12 read file wrms name=(rms.data) unit=14 wovr #crd=100 #ste=1000 step=0.002 action below is the standard output fluctuation> ~ fluctuation> ~ DYNAMICS TEST fluctuation> ~ fluctuation> file conn name=(wcon.mbco) unit=10 read fluctuation> file rcrd name=(MBCO.MIN.CRD) unit=11 read getcrd> * title for CHARMM coordinates getcrd> * fluctuation> file rdyc name=(MBCO.DCD) bina unit=12 read fluctuation> file wrms name=(rms.data) unit=14 wovr fluctuation> #crd=100 #ste=1000 step=0.002 fluctuation> action READING DYNAMICS CRD, STRUCTURE # 1 rotation matrix .9999999999927429 2.642858385206576E-06 2.744106114227484E-06 -2.642869632723640E-06 .9999999999881071 4.098891828923495E-06 -2.744095281424080E-06 -4.098899081242629E-06 .9999999999878348 After overlap rms is = .1800677569844866 READING DYNAMICS CRD, STRUCTURE # 2 rotation matrix .9999999997731057 1.219989243025524E-05 1.746289086074909E-05 -1.220010823773631E-05 .9999999998492188 1.235800795341013E-05 -1.746274009178233E-05 -1.235822099971166E-05 .9999999997711643 After overlap rms is = .3206856287223989 READING DYNAMICS CRD, STRUCTURE # 3 rotation matrix .9999999996325165 4.004725312505257E-06 2.681284716665732E-05 -4.004024638479703E-06 .9999999996505422 -2.613203795784942E-05 -2.681295180890358E-05 2.613193058893445E-05 .9999999992990935 After overlap rms is = .4447835002218578 READING DYNAMICS CRD, STRUCTURE # 4 rotation matrix .9999999994612092 3.091141573821965E-05 -1.104836711908474E-05 -3.091149265264087E-05 .9999999994980089 -6.961487758115787E-06 1.104815192408750E-05 6.961829275893387E-06 .9999999999147361 After overlap rms is = .5234805210556039 READING DYNAMICS CRD, STRUCTURE # 5 rotation matrix .9999999966409292 7.807992912594330E-05 -2.493321283645113E-05 -7.807987782327046E-05 .9999999969496477 2.058565041526466E-06 2.493337349297047E-05 -2.056618252384680E-06 .9999999996870484 After overlap rms is = .5859792965449386 READING DYNAMICS CRD, STRUCTURE # 6 rotation matrix .9999999951942369 7.925194567436919E-05 5.771182729550141E-05 -7.925096753584148E-05 .9999999967159779 -1.695075879757970E-05 -5.771317048658755E-05 1.694618499797928E-05 .9999999981910089 After overlap rms is = .6540930897502345 READING DYNAMICS CRD, STRUCTURE # 7 rotation matrix .999999985564983 1.070426953304924E-04 1.319541357243832E-04 -1.070417189140322E-04 .9999999942436046 -7.406707871853801E-06 -1.319549277988066E-04 7.392583167415887E-06 .9999999912666236 After overlap rms is = .7424947261572076 READING DYNAMICS CRD, STRUCTURE # 8 rotation matrix .9999999929847465 1.180966651069081E-04 -9.147913594558659E-06 -1.180951365161499E-04 .9999999790959805 1.669178782365706E-04 9.167625848065242E-06 -1.669167967414957E-04 .9999999860273685 After overlap rms is = .8007353643772485 READING DYNAMICS CRD, STRUCTURE # 9 rotation matrix .9999999957592453 7.789672531365662E-05 -4.912850703480894E-05 -7.787871896783624E-05 .9999999298390524 3.664106927434027E-04 4.915704578100505E-05 -3.664068651243663E-04 .9999999316647944 After overlap rms is = .8447578167717348 READING DYNAMICS CRD, STRUCTURE # 10 rotation matrix .9999999994376918 3.231737448911736E-05 -8.955677125409290E-06 -3.231343068263637E-05 .9999999026690097 4.400202425462557E-04 8.969896552749739E-06 -4.400199529101490E-04 .9999999031509859 After overlap rms is = .8771241859329588 CONN.FIL This file defines the connectivity of a specific molecule. It is a formatted output of the program con that can be used in energy calculations, dynamics, etc. knowledge on this file is not essential for normal use of the program, however it is useful to know what it includes for debugging purposes and general understanding of how the program works. Below there is a real connectivity file generated for alanine dipeptide Clarifying comments (that are not included in the original file) are given between {*** and ***}. Also added empty lines separating comments and real stuff to increase readability. ~ CONNECTIVITY FILE FOR MOLECULES: {*** totmon is the toal number of monomers in the generated macromolecule npt is the total number of particles nb is the total number of bonds nangl is the total number of angles ntors is the total number of torsions nimp is the number of improper torsions totex total number of exclusions (no vdw and electrostatic for 1-2 1-3 interactions totspe is the total number of special 1-4 interactions lestyp is the number of different les type particles available NBULK is the number of BULK (macromolecules generated here) ***} ~ totmon npt nb nangl ntors nimp totex totspe lestyp NBULK 4 15 14 21 9 5 35 22 0 1 {*** Below the name(s) of the BULK are given ***} ALA1 {*** That pointer(s) are the last particle(s) of the macromolecule(s). This example has only one and this is also the last atom. ***} ~ Pointers to last particle of BULK 15 {*** The name NTER is a special monomer for N terminal, CTER is the caroxyl terminus monomer. Hence the molecule below is dialanine ***} ~ Monomer names NTER ALA ALA CTER {*** Given the last atom of all the monomers ***} ~ Pointers to last particle of monomer 2 8 14 15 {*** List of particle properties: pt - a running index over the particle numbers mono - a running index over the monomer index. ptid - index of particle type, can be traced back to the PTNM of the property file. lesid - index of LES particles. 0-normal particles ; 1-... different LES groups. Particles of the same group do not see each other ptnm - particle unique name (withina monomer) ptms - particle mass (g/mol units) ptchg - particle charge (charge of electron units) epsgm6 & epsgm12 - for vdw interactions 4*epsilon*sigma^n (n=6,12) ptwei - weight of LES particle (between zero and one) ***} ~ Properties of particles list : ~pt mono ptid lesid ptnm ptms ptchg epsgm6 epsgm12 ptwei 1 1 2 0 HX2 1.00 0.33000 0.00000E+00 0.00000E+00 0.10000E+01 2 1 2 0 HX3 1.00 0.33000 0.00000E+00 0.00000E+00 0.10000E+01 3 2 1 0 N 14.00 -0.30000 0.28308E+02 0.97175E+03 0.10000E+01 4 2 2 0 H 1.00 0.33000 0.00000E+00 0.00000E+00 0.10000E+01 5 2 3 0 CA 13.00 0.31000 0.37698E+02 0.20686E+04 0.10000E+01 6 2 12 0 CB 15.00 0.00000 0.47821E+02 0.28586E+04 0.10000E+01 7 2 9 0 C 12.00 0.50000 0.34176E+02 0.18022E+04 0.10000E+01 8 2 10 0 O 16.00 -0.50000 0.23769E+02 0.61644E+03 0.10000E+01 9 3 7 0 N 14.00 -0.57000 0.22247E+02 0.76370E+03 0.10000E+01 10 3 8 0 H 1.00 0.37000 0.00000E+00 0.00000E+00 0.10000E+01 11 3 6 0 CA 13.00 0.10000 0.19961E-02 0.10220E-05 0.10000E+01 12 3 12 0 CB 15.00 0.00000 0.47821E+02 0.28586E+04 0.10000E+01 13 3 4 0 C 12.00 0.70000 0.34176E+02 0.18022E+04 0.10000E+01 14 3 5 0 O 16.00 -0.80000 0.23769E+02 0.61644E+03 0.10000E+01 15 4 5 0 OX2 16.00 -0.80000 0.23769E+02 0.61644E+03 0.10000E+01 {*** The indices ib1 and ib2 refer to the pt running index, kbond is in kcal/mol angstrom^-2 req in angstrom ***} ~ Bonds list: ~ ib1 ib2 kbond req 1 3 434.0000 1.0100 2 3 434.0000 1.0100 3 5 337.0000 1.4490 3 4 434.0000 1.0100 5 6 260.0000 1.5260 5 7 317.0000 1.5220 7 9 490.0000 1.3350 7 8 570.0000 1.2290 9 10 434.0000 1.0100 9 11 337.0000 1.4490 11 12 260.0000 1.5260 11 13 317.0000 1.5220 13 15 656.0000 1.2500 13 14 656.0000 1.2500 {*** similar to bonds. Units: kangle - kcal/mol radians^-2 ; angleq - degrees ***} ~ Angles list: ~ iangl1 iangl2 iangl3 kangl angleq 1 3 2 35.00000 109.50000 1 3 4 35.00000 109.50000 1 3 5 35.00000 109.50000 2 3 4 35.00000 109.50000 2 3 5 35.00000 109.50000 4 3 5 35.00000 109.50000 3 5 6 35.00000 109.50000 3 5 7 63.00000 110.10000 6 5 7 63.00000 111.10000 5 7 8 80.00000 120.40000 5 7 9 70.00000 116.60000 8 7 9 80.00000 122.90000 7 9 10 35.00000 119.80000 7 9 11 50.00000 121.90000 10 9 11 38.00000 118.40000 9 11 12 80.00000 109.50000 9 11 13 63.00000 110.10000 12 11 13 63.00000 111.10000 11 13 14 65.00000 117.00000 11 13 15 65.00000 117.00000 14 13 15 80.00000 126.00000 {*** as bonds and angles. See prop.file for parameter explan ***} ~ Torsions list: ~ itor1 itor2 itor3 itor4 period ktors1 ktors2 ktors3 phase 3 5 7 8 3 0.0000 0.0000 0.1000 -1.000 5 7 9 10 2 0.0000 2.5000 0.0000 -1.000 5 7 9 11 2 0.0000 2.5000 0.0000 -1.000 8 7 9 10 2 0.0000 2.5000 0.0000 -1.000 8 7 9 11 2 0.0000 2.5000 0.0000 -1.000 9 11 13 14 3 0.0000 0.0000 1.0000 1.000 9 11 13 15 3 0.0000 0.0000 1.0000 1.000 12 11 13 14 3 0.0000 0.0000 1.0000 1.000 12 11 13 15 3 0.0000 0.0000 1.0000 1.000 {*** As above ***} ~ Improper torsion properties: ~iimp1 iimp2 iimp3 iimp4 kimp impeq 5 3 7 6 .5500000000E+04 .3525999832E+02 7 5 8 9 .1000000000E+05 .0000000000E+00 9 11 7 10 .4500000000E+04 .0000000000E+00 11 9 13 12 .5500000000E+04 .3525999832E+02 13 11 14 15 .1000000000E+05 .0000000000E+00 {*** organized as follows: for particle i there are n other particles that do not interact with it because they are one bond (1-2) or two bonds (1-3) apart. they are j1,j2,....,jn. Everything is written as i n j1 j2 j3 ....jn ***} ~ Exclusion list 1-2 1-3, set as followed: ~ atom number, number of exclusions and list 1 4 3 2 4 5 2 3 3 4 5 3 4 5 4 6 7 4 1 5 5 4 6 7 8 9 6 1 7 7 4 9 8 10 11 8 1 9 9 4 10 11 12 13 10 1 11 11 4 12 13 14 15 12 1 13 13 2 15 14 14 1 15 {*** organized in the same way as the exclusion list ***} ~ Special list 1-4 set as followed: ~ atom number, number of exclusions and list 1 2 6 7 2 2 6 7 3 2 8 9 4 2 6 7 5 2 10 11 6 2 8 9 7 2 12 13 8 2 10 11 9 2 14 15 10 2 12 13 12 2 14 15 {*** This is it! ***} MONO.FIL The monomer file is where the connectivity of the particles in a monomer is listed and where the rules for joining monomers are defined. It is an input to the "con" program in order to generate the connectivity file for the complete molecule. The structure of the file is as follows: The top of the file must be (or a comment line that have ~ ) MONO LIST after the title the different monomers are listed. Each monomer starts with the line MONO=(NAME) #prt=5 chrg=0. where NAME is the name of the monomer type that can be at most four characters (all character assignments must be closed in brackets. This include of course also the monomer name). #prt is the number of particles in the monomer. This includes also virtual particles used to link the monomer to next or previous monomers.. chrg is the charge of the total monomer INCLUDING the virtual particles. The virtual particles are included since their type is the one that will be finally used (see below). The total charge is used only for test purposes. After that line ca list of unique names of particles (to that monomer) and their types is provided: ~ unique particle to monomer particle type connectivity action UNIQ=(UNAM) PRTC=(PTYP) HERE UNIQ=(B) PRTC=(BTYP) UNIQ=(C) PRTC=(CTYP) NEXT where the UNIQ command assigns a unique particle name (unique to that monomer). The name can be at most four characters. The name must appear only once for a given type of connectivity action (see below), otherwise bond assignments will be ambiguous. The PRTC defined the particle type and it is matched against the PNAM data from the property file. Failure to match particle types results in program termination. The connectivity action provides the necessary data to join monomers. The following keywords are available for the connectivity action (note that only four characters from a keyword are actually used): HERE (or blank) - This is a normal particle of the present monomer NEXT - This particle belongs to the next monomer, when the Molecule will be Formed (MF), do the following: Change the Type of the Particle (CTP) in the next monomer to the one defined in the present monomer, Transfer all Bonds (TB) in the present monomer to the one that it is connected to. PREVIOUS - The particle belongs to the previous monomer, when MF do CTP and TB DNXT - At MF delete particle and all its bond in next monomer. DPRV - At MF delete particle and all its bond in previous monome. The unique particle list ends (as usual) by DONE after the particle list, a bond list is provided. A Bond is given by two unique particle names and a dash in between (e.g. A-B is a bond between A and B). Note that special particles for which the connectivity action is different from HERE are denoted by *. Example following the definition of particles above is. BOND B-C* UNAM-C* B-UNAM DONE Special particles must come second in bond and a monomer cannot be used to define a bond between two special particles, one of the particles must be HERE. Note that the * is really needed to avoid ambiguity. Since it is possible to have (for example), particle A as HERE and also particle A as PREV. This also means that there are some restrictions on the connectivity. Currently it is not possible to make a reference from a given monomer to the same particle name at PREVious and NEXT monomers. It is hard to imaging however a case in which it is truly needed. Note also that cases in which a NEXT particle is defined but not found are possible, a warning will be issued but it is not fatal. In fact it can be quite convenient to define peptide link as the carbonyl carbon attached to the NEXT nitrogen. This however fails at the C terminus. The program is set in this way. The virtual nitrogen at the C terminus is deleted which bring everything back to normal, a warning on that nitrogen is however issued, that warning should be ignored. The angles, torsions improper torsions are generated once the complete molecule is formed and they are generated in a complete way. I.e. all possible angles, torsions and improper torsions are formed. Some torsions are then eliminated (with zero energy contribution). One consequence is that the possible bonds angles and improper torsions MUST be defined in the property file. If torsion is not found a yellow alert is issued (non-fatal warning) and that torsion is ignored. The file ends in the traditional way, i.e. *EOD Finally we give a complete example to define alanine MONO=(ALA) #prt=7 chrg=-0.57 UNIQ=(N) PRTC=(NH) UNIQ=(H) PRTC=(HN) UNIQ=(CA) PRTC=(CAH) UNIQ=(CB) PRTC=(CH3) UNIQ=(C) PRTC=(CO) UNIQ=(O) PRTC=(OC) UNIQ=(N) PRTC=(NH) NEXT DONE BOND C-O C-N* C-CA CA-CB CA-N N-H DONE An example for a monomer file can be found /store/us/ron/moil.mop/ALL.MONO POLY.FIL This (short) file gives a list of monomer for a given molecule (sequence). It is used in the program con. Format is simple (I think). MOLC is the molecule name (called BULK in the conn file) four characters at most. #mon is the number of monomers. The lines after, includes the monomer names. The number of monomers found should match the number of monomers declared (#mon). MOLC=(BIG) #mon=3 NTER ALA CTER *EOD To define LES monomers "mult" should be used. For example, to link NTER ALA ALA(les with ten copies) ALA CTER we use MOLC=(BIG) #mon=14 NTER ALA mult ALA 10 ALA CTER *EOD when a "mult" is found in a sequence, the next expression is the monomer name, that monomer will be multiplied the number that comes immediately after. This means that the connectivity file will have all the appropriate id's and normalization (if all the bugs were eliminated...). Currently it is not possible to multiply only a part of monomer (unless it is defined as a monomer!). The number of possible monomers is unlimited though one should be careful not to exceed to maximum number of bonds to a particle with is defined in LENGTH.BLOCK Another option added!!. You may specify the index of your LES monomer. This means that also two different monomers may be set NOT to see each other (for example two different sidechains connected to the same CA carbon - Bob Goldstein idea) The polymer file MOLC=(BIG) #mon=13 NTER ALA mult ALA 5 mult ALA 5 ltyp 6 CTER *EOD will set the lesid of the second mutiplied residue to 6 (instead of the usual increment by 1 that will set it to two. "ltyp 0" will set the lesid to be the same as previous monomer. PROP.FIL The property file is where the parameters of the particles, bonds, angles, torsions and improper torsions are kept. It is used as an input for the connectivity program which build the molecular connectivity file The property file is build from sequential sections which MUST come in the following order: PRTC - individual particle properties BOND - bond parameters ANGLE - angle parameters TORSION - torsion parameters IMPROPER - improper torsion parameters It is possible not to provide all the information i.e. a property file with PRTC only is legal, however PRTC and ANGLE is not. If you provide only PRTC the program will issue a warning (yellow alert), ignore it unless you want to provide bonds and something happened. Each subsection (e.g. PRTC, BOND, TORSION) must end with DONE The file must end with *EOD The DONE and *EOD are general termination features used in other datafile. Another general feature shared between different data files is the comment line. ~ ANYWHERE in the line makes it a comment. This line is ignored by the program. Below details on the syntax are provided: The explanations will be written as comment lines as in a "real" property file ~ This is a first line of a property file. The first exe line must be PRTC ~ The following line lists properties of an individual particle ~ name mass charge epsilon sigma PNAM=(NX) PMAS=14. PCHG=-0.3 PEPS=0.170 PSGM=3.250 ~ The example above provides the data for the particle type NX (Nitrogen ~ of the N-terminal. characters (like PNAM - the name of particle ~ type) must be enclosed in brackets. Each of the expressions (i.e. ~ A=B) must be seprated from other expressions by space(s). No spaces ~ within an expression are allowed. ~ espilon and sigma are the van der Waals well depth ~ and the hard core radius respectively. Obvioulsy this type of ~ line is repeated as needed for different particle types. ~ The data base for particle properties is based on the OPLS potential ~ Jorgensen and Tirado-Rives JACS 110,1657(1988) ~ ~ Now end the particle part by DONE DONE ~ ~ The next part lists the bond properties. The first Bond line must be BOND ~ Bond energy is set to be K(r - req)^2 ~ Below we provide the names of the two particle types, the force constant ~ (in kcal/mol angstrom^-2) and the equilibrium distance in angstrom ~ Note the different style of i/o different expression are still ~ separated by spaces but no equality is used. This requires the data ~ to be placed in exactly the same order. I.e. do not exchange equilibrium ~ position and force constant. ~ The covalent part of the potential (excluding improper torsions) ~ is taken from AMBER ~ Weiner et al JACS 106,765(1984) ~ particle particle force-constant equilbrium distance NX HX 434.0 1.01 CANX NX 337.0 1.449 ~ ~ Pictorialy NX-HX ~ end the BOND with DONE DONE ~ ~ Angles are similar to bonds in format style ANGLE ~ K (theta -theta(eq))^2 ~ name name name K(kcal/mol radians^-2) theta(eq) (degrees) HX NX HX 35.0 109.5 ~ ~ Pictorially HX-NX-HX DONE ~ ~ And here are the torsions. The format style is similar to BOND and ANGLE TORSION ~ (Pictorially CAH-CO-NH-CAH) ~ however the energy function is more complex: ~ E = sum k(n)*(1 + cos(n*phi+gamma) ~ (gamma should be a function of n too and will be added to the program ~ soon). Currently the format is ~ name name name name k(1) k(2) k(3) n cos(gamma) CAH CO NH CAH 0.0 2.5 0.0 2 -1.0 ~ There is an option in TORSION (only) to use a wild card by X, e.g. X CANX CX X 0.0 0.0 0.0 3 0.0 ~ where X means "any atom". ~ **** ALL TORSIONS MUST BE DEFINED IN THE PROPERTY FILE *** ~ However in many cases the energy is set identically to zero. ~ This is done by setting cos(gamma)=0. When the program matched ~ this torsion, it is skipped and NOT included finally in the ~ connectivity file ~ DONE ~ ~ Improper torsions are four body interactions in which one atom is sitting ~ in the center, Pictorially ~ B ~ | ~ A ~ / \ ~ C D IMPROPER ~ The internal degree of freeom - phi, is the angle between the normal ~ to the ABC plane and the normal to the BCD plane. To obtain consistent ~ values A must be first and D must be last ~ ~ The energy function is rather messy... ~ If the equilibrium angle is far from zero we use simply harmonic term ~ E = K1(phi-phi(eq))^2 (K1 kcal/mol radian^-2 ; phi(eq) degrees) ~ If the equilibrium angle equals zeo then the above energy expression ~ is singular, we therefore use ~ E = K2(cos(phi) - cos(phi(eq))^2 (K2 kcal/mol ; phi(eq) degrees) ~ Note that the units for K is different in both cases. Note also ~ that at phi(eq)=0, Taylor exapnsion shows that E is quartic in phi ~ therfore to maintain comparable restoring force K2 > k1 ~ ~ The atom in the center is always first, the last atom must also be chosen ~ with care since it determines the sign. The other two in the middle ~ can be interchange ~ name name name name K1/K2 phi(eq) CANX NX CO CH3 55.0 35.26 DONE ~ End it all *EOD An example for a property file can be found in the directory bully:/store/us/ron/moil.mop/ALL.PROP SUB.LIST connect.f: program con doing the connectivity chain/chain.f: program chain doing chain simulation dynamics/dyna.f: program dynamics doing dynamics minimize/mini.f: program mini doing energy minimization pot/energy.f: program pote calculating current POTential Energy addbond.f: subroutine addbond(ubond) add bonds to existing BULK angtor.f: subroutine angtor() make angles torson and improper torsions angtor.f: subroutine pickup2(n, atom, arr1, arr2, out) angtor.f: subroutine sort1(n, ra) angtor.f: subroutine sort2(n, ra, rb) angtor.f: subroutine sort4(n, ra, rb, rc, rd) angtor.f: subroutine piksrt(n, arr) dimer.f: subroutine dimer(ityp1,ityp2,mono_pt,updown,nprt,nbnd,iden, joining two monomers:bonds, etc fillpt.f: subroutine fillpt(ipart1,ityp1,ptnm_m1,id_m1,ib1_m1,ib2_m1 fill particles in BULK getcrd.f: subroutine getcrd(ucrd,type) read coordinates getmoname.f: subroutine getmoname(name,imult,less,empty,copies read monomer names from poly file lego.f: subroutine lego(upoly,mono,nprt,mono_pt,nbnd, top routine for ploymerization process monomers.f: subroutine monomers(ucon,mono,nprt,mono_pt,nbnd,basenm, read and organize mono file overlap.f: subroutine overlap(x,y,z,x1,y1,z1) overlap two structures for minimum mass weighted rms pick.f: subroutine pick(ipick,ngrp) pick.f: subroutine pick2(ipick,ngrp) pick atoms (not working yet) propert.f: subroutine prprt(uprm,basenm,basems,basechg, read and organize property file putcrd.f: subroutine putcrd(ucrd,type) write coordinates (CHARMm format formatted) rconn.f: subroutine rconn(urcon) read connectivity file search.f: subroutine srch_imp(indx1,indx2,indx3,indx4,i1,i2,i3,i4, search for improper torsion in data_base vcomprs.f: subroutine cmprs(vector,away,length) compress one vector vcomprs.f: subroutine cmprs2(vector1,vector2,away,length) compress two vectors vcopy.f: subroutine vdcopy(v1,v2,length) copy double precision vector vcopy.f: subroutine vicopy(v1,v2,length) copy an integer vector vcopy.f: subroutine vccopy(v1,v2,length,size) copy a character vector wconn.f: subroutine wconn(uwcon) write connectivity file chain/chndyn.f: subroutine chndyn(temp,tmpr,step,tfac,igrid,npri,nstep,nsvel, do CHAIN dynamics chain/comc.f: subroutine comc(r0,dmass,natom,pointr, generate rigid body constraints chain/crbm.f: subroutine crbm(r,scalar,sigma,grdcmx,grdcmy,grdcmz, project out rigid body motion chain/orthg.f: subroutine orthg(X,Y,n,amovr,sigma,n1,n2) orthogonalize vectors with weight chain/sds.f: subroutine sds(sener,dv,r,d0,e0,e1,bnbnd,bimag,nselec,pointr, calculate derivative and energy of the whole chain dynamics/getvel.f: subroutine getvel(uvel,type) get velocity from coordinate file dynamics/info.f: subroutine info(unit,jstep,nofreez,inofrz) write info on dynamics progress dynamics/picktemp.f: subroutine picktemp pick multiple temperature (NOT WORKING YET) dynamics/rdyncrd.f: subroutine rdyncrd(urcrd,ncoor,inofrz,nofreez) read dynamics coordinate dynamics/shake.f: subroutine shake(nc, ishak1, ishak2, dissq, shake bonds (Not implemented yet) dynamics/shakept.f: subroutine shakept() dummy shake routine (not working yet) dynamics/shakinit.f: subroutine shakinit(nofreez,inofrz) initialize shake database (Not implemented yet) dynamics/velinit.f: subroutine velinit(irand,tempi,ntemp,tpo) initialize velocities from Boltzmann distribution dynamics/wdyncrd.f: subroutine wdyncrd(uwcrd,nfrm,ncoor,inofrz,nofreez) write CHARMm dynamics coordinate file dynamics/wdynvel.f: subroutine wdynvel(uwvel,nfrm,nvelo,inofrz,nofreez) write CHARMm dynamics velocity file house.f: subroutine house householder diagonalization routine intrprtr/alert.f: subroutine alert(func,length1,message,length2,level) Error messenger intrprtr/brkpair.f: subroutine brkpair(exp1,l1,exp2,l2,special) breaking pairs (used to read bonds) intrprtr/echo.f: subroutine echo(sbr,sbrl,line,linel) echoing command line intrprtr/get4c.f: subroutine get4c(c4,empty) get character*4 from command line intrprtr/rline.f: subroutine rline(name,namel,unit) read current command line minimize/eforce0.f: subroutine eforce0(idim,xx,etotal,g) initiate energy call from the minimizer ******THE PATH ROUTINES ARE NOT WORKING YET **** path/pat17.f: subroutine PAT17(COMLYN,COMLEN,DS,ATOMPR,RW,EN,DC,WORK pot/ebond.f: subroutine ebond() bond energy pot/eforce.f: subroutine eforce(cutoff) top of energy routines pot/eimphi.f: subroutine eimphi() improper energy pot/enonb.f: subroutine enonb() van der Waals and electrostatic interactions pot/ephi.f: subroutine etors() torsion energy pot/etheta.f: subroutine etheta() angle energy pot/exlist.f: subroutine exlist(l,m,flag) sort and exclude those particles in exclusion list pot/exlist.f: subroutine incl14(l,m,flag14) 1-4 list handler pot/init_wre.f: subroutine init_wre(uwene,cutoff) initiate energy output pot/nbond.f: subroutine nbond(pointm,nblistm) create atom neighbour list pot/nbondm.f: subroutine nbondm(cutoff) create monomer-monomer neighbour list pot/wener.f: subroutine wener(uwene) write energy to uwene file of.f: function of() Open File search.f: integer function search(indx1,indx2,indx3,indx4, search and match abgles or torsions against data-base chain/hotchain.f: function hotchaf(vx,vy,vz,ptms,inofrz) provide current temeprature of chain dynamics/hot.f: function hot(vx,vy,vz,ptms,nofreez,inofrz) provide current temeprature of dynamics dynamics/velinit.f: function ran1(irand) random number generator (poor one) dynamics/velinit.f: function ran2(irand) high quality random number generator intrprtr/find.f: logical function find(expr) if find expr in a line find=.true. intrprtr/fopen.f: logical function fopen(unit) if unit is opened fopen=.true. intrprtr/getc.f: character*(*) function getc(expr,def,lcpext) get a character with unspecified length intrprtr/geti.f: integer function geti(expr,def) get integer from expression in line intrprtr/getr.f: function getr(expr,def) get real number from expression in line intrprtr/intgr.f: function intgr() get the integer from next expression intrprtr/number.f: function number() get real number from next expression TUTORIAL.SHORT WELCOME TO MOIL! (A tutorial by Danuta Rojewska) We'll begin by building a benzene molecule "by hand", that is by setting up and using your own data base. The purpose of this exercise is to show you that it is a relatively easy task to "invent" and generate molecules other than amino-acids on your own. A MOIL "molecule" consists of one or more segments (or residues) that we call monomers, alligned in a specified sequence. In order to generate one you need to provide MOIL with 3 data files which contain: 1. the PROPerties of the different particle types which make up a monomer ( mass, charge, etc.) and the covalent structure of the monomer (force constants and equilibruim values for bonds, angles and torsions). NOTE: The particle type parameters come from the OPLS paper by Jorgensen and Tirado-Rives (JACS 110,1657(1988)), The covalent parameters come from AMBER paper Weiner et al. JACS 106,765(1984) (The names of the particles are different). 2. the description of the individual MONOmers which you may link together, if you want, in the next step, 3. the POLYmerization info, where you tell the program what sequence of monomers you want. You also need to tell MOIL where to write the output connectivity file (WCON) if you want to use this molecule later. Here is the input file which generates the connectivity for benzene. The lines starting with ~ are comment lines, and won't be executed. ~ ~ input file ~ ~ tell MOIL where to find your files and what to do with them: file prop name=(BENZ.PROP) unit=10 read file mono name=(BENZ.MONO) unit=11 read file poly name=(BENZ.POLY) unit=15 read file wcon name=(BENZ.WCON) unit=13 wovr ~ now prod the program to do some more serious activity: action *EOD ~ ~ "wovr" means write over an existing file if there is one. Other option is "write". This one will bounce you if you try to write on a file with the same name. Here is what you need to put in those files: BENZ.PROP file ------------------------------------------------------------------------- ~ first, the info about the individual particle types: PRTC ~ here C and H are treated as a united atom: PNAM=(CH) PMAS=13.0 PCHG=0.00 PEPS=0.110 PSGM=3.750 ~ ~ end this subsection: DONE ~ BOND ~ name name force constant equilibrium value CH CH 469.0 1.40 DONE ~ ANGLE CH CH CH 85.0 120.0 DONE ~ TORSION ~ If you are not sure what a torsion angle is: ~ For a set of 4 atoms connected in a sequence, ABCD, ~ a proper torsion is the angle between the planes of ABC and BCD ~ B D ~ / \ / ~ A C ~ ~ NOTE: The torsion parameters are derived from AMBER ~ but MOIL uses a different notation: ~ E = SUM V(n)/2( 1 + cos(n*phi+gamma) ~ as k(n)=V(n)/2 and phase=cos(gamma) ~ So the MOIL data structure is: ~ name name name name k(1) k(2) k(3) n cos(gamma) ~ CH CH CH CH 0.0000 5.3000 0.0000 2 -1.0 ~ ~ the k's are the Fourier coefficients V(n)/2 ~ and -1 in this case corresponds to angle ~ gamma equal to 180 degrees. ~ DONE ~ IMPROPER ~ Now if you are not sure what an improper torsion is, here is the ~ atomic arrangement it describes. Basically, it is the angle between ~ planes ACD and BCD, giving one an idea as to how much atom C goes out ~ of the plane defined by atoms ABD. ~ ~ A B ~ \ / ~ C ~ | ~ D ~ Thus, benzene doesn't really have improper torsions in our rendition ~ (i.e., for CH taken as a unified atom) ~ ~ taken from CHARMM paper (Brooks et al. J. Comput. Chem ...) CH CH CH CH 25.0 0.0 DONE ~ ~ end of data for this section: *EOD ---------------------------------------------------------------------- BENZ.MONO file: ---------------------------------------------------------------------- MONO LIST ~ we will only define one monomer, but more can be listed ~ (start with MONO= finish with DONE, just like for this one); ~ ~ this monomer has 6 particles and zero total charge: MONO=(BENZ) #prt=6 chrg=0. ~ ~ the unique particle names differentiate the atoms within the ~ monomer, even though they may be of the same type: UNIQ=(C1) PRTC=(CH) UNIQ=(C2) PRTC=(CH) UNIQ=(C3) PRTC=(CH) UNIQ=(C4) PRTC=(CH) UNIQ=(C5) PRTC=(CH) UNIQ=(C6) PRTC=(CH) DONE ~ don't forget to state which particles are linked together: BOND C1-C2 C2-C3 C3-C4 C4-C5 C5-C6 C6-C1 DONE ~ ~ tell the program that's it for this file: ~ (the angles, torsions, and improper torsions are generated automatically, ~ so there is no need to specify them in the monomer file) *EOD ~ BENZ.POLY file ------------------------------------------------------------------------ ~ name your molecule and tell MOIL the number of monomers you ~ want linked: MOLC=(BENZ) #mon=1 ~ ~ and in what sequence: BENZ *EOD --------------------------------------------------------------------- To run the input: Set your path variable to include moil.src.001/exe and type (assuming your file is called "input"): con < input > output This will activate the part of MOIL that deals with connectivity. The output should give you the following: -------------------------------------------------------------------- CONN> file prop name=(BENZ.PROP) unit=10 read CONN> file mono name=(BENZ.MONO) unit=11 read CONN> file poly name=(BENZ.POLY) unit=15 read CONN> file wcon name=(BENZ.WCON) unit=13 wovr CONN> action propert> PRTC propert> PNAM=(CH) PMAS=13.0 PCHG=0.00 PEPS=0.110 PSGM=3.750 propert> DONE *** DONE WITH PARTICLE PROPERTIES *** propert> BOND propert> CH CH 469.0 1.40 propert> DONE *** DONE WITH BOND PROPERTIES *** propert> ANGLE propert> CH CH CH 85.0 120.0 propert> DONE *** DONE WITH ANGLE PROPERTIES *** propert> TORSION propert> CH CH CH CH 0.0000 5.3000 0.0000 2 -1.0 propert> DONE *** DONE WITH TORSION PROPERTIES *** propert> IMPROPER propert> CH CH CH CH 25.0 0.0 propert> DONE *** DONE WITH IMPROPER TORSIONS *** *** Number of parameters read: Unique particles, bonds, angles, torsions, improper torsions 1 1 1 1 1 monomers> MONO LIST monomers> MONO=(BENZ) #prt=6 chrg=0. monomers> UNIQ=(C1) PRTC=(CH) monomers> UNIQ=(C2) PRTC=(CH) monomers> UNIQ=(C3) PRTC=(CH) monomers> UNIQ=(C4) PRTC=(CH) monomers> UNIQ=(C5) PRTC=(CH) monomers> UNIQ=(C6) PRTC=(CH) monomers> DONE monomers> ~ monomers> BOND monomers> C1-C2 C2-C3 C3-C4 C4-C5 C5-C6 C6-C1 monomers> DONE monomers> ~ monomers> *EOD LEGO> MOLC=(BENZ) #mon=1 LEGO> BENZ *** MOLECULE BENZ ASSEMBLED *** LEGO> *EOD ---------------------------------------------------------------------- Now, just in case, look at the resulting connectivity file, to see if it looks reasonable. Here is what my BENZ.WCON looks like: ---------------------------------------------------------------------- ~ CONNECTIVITY FILE FOR MOLECULES: ~ totmon npt nb nangl ntors nimp totex totspe lestyp NBULK 1 6 6 6 6 0 15 0 0 1 BENZ ~ Pointers to last particle of BULK 6 ~ Monomer names BENZ ~ Pointers to last particle of monomer 6 ~ Properties of particles list : ~pt mono ptid lesid ptnm ptms ptchg epsgm6 epsgm12 ptwei 1 1 1 0 C1 13.00 .00000 .34980E+02 .18446E+04 .10000E+01 2 1 1 0 C2 13.00 .00000 .34980E+02 .18446E+04 .10000E+01 3 1 1 0 C3 13.00 .00000 .34980E+02 .18446E+04 .10000E+01 4 1 1 0 C4 13.00 .00000 .34980E+02 .18446E+04 .10000E+01 5 1 1 0 C5 13.00 .00000 .34980E+02 .18446E+04 .10000E+01 6 1 1 0 C6 13.00 .00000 .34980E+02 .18446E+04 .10000E+01 ~ Bonds list: ~ ib1 ib2 kbond req 1 2 469.0000 1.4000 1 6 469.0000 1.4000 2 3 469.0000 1.4000 3 4 469.0000 1.4000 4 5 469.0000 1.4000 5 6 469.0000 1.4000 ~ Angles list: ~ iangl1 iangl2 iangl3 kangl angleq 2 1 6 85.00000 120.00000 1 2 3 85.00000 120.00000 2 3 4 85.00000 120.00000 3 4 5 85.00000 120.00000 4 5 6 85.00000 120.00000 1 6 5 85.00000 120.00000 ~ Torsions list: ~ itor1 itor2 itor3 itor4 period ktors1 ktors2 ktors3 phase 3 2 1 6 2 .0000 5.3000 .0000 -1.000 2 1 6 5 2 .0000 5.3000 .0000 -1.000 1 2 3 4 2 .0000 5.3000 .0000 -1.000 2 3 4 5 2 .0000 5.3000 .0000 -1.000 3 4 5 6 2 .0000 5.3000 .0000 -1.000 1 6 5 4 2 .0000 5.3000 .0000 -1.000 ~ Exclusion list 1-2 1-3, set as followed: ~ atom number, number of exclusions and list 1 5 2 6 3 5 4 2 4 3 6 4 5 3 3 4 5 6 4 2 5 6 5 1 6 --------------------------------------------------------------------------- The description of all the variables can be found in the file conn.file Now we can check if our molecule is comfortable and help it relax: Use input like the one below for the energy calculation. All you need here is 1. the connectivity file (*.WCON), and 2. the atomic coordinates in the CHARMM format, (*.CRD ) A detailed description of the setup options for the energy calculation can be found in the file energy.prog For a simple molecule like this we can turn off all energy terms but the bond and angle energies (energy.prog will give you the full list of commands) but they will be required for peptides and proteins. Here is my coordinate file for benzene (BENZ.CRD): can be standard CHARMM Coordinate file Organized as follows Title: * text Title: * total number of particles Crd format: particle no. mono.no. mono.name atom.name x y z junk junk fourth --------------------------------------------------------------------------- * benzene (CHARMM format) * 6 1 1 BENZ C1 10.10440 46.18651 12.13924 BENZ 1 13.01900 2 1 BENZ C2 11.12320 46.28111 11.21186 BENZ 1 13.01900 3 1 BENZ C3 10.96320 47.07198 10.09118 BENZ 1 13.01900 4 1 BENZ C4 9.78651 47.76534 9.89150 BENZ 1 13.01900 5 1 BENZ C5 8.77222 47.68540 10.82364 BENZ 1 13.01900 6 1 BENZ C6 8.93509 46.90089 11.95008 BENZ 1 13.01900 ---------------------------------------------------------------------------- ~ ~ input file: ~ ~ energy calculation: ~ file rcon name=(BENZ.WCON) unit=10 read file rcrd name=(BENZ.CRD) unit=11 read file wene name=(BENZ.ENE) unit=13 wovr rmax=9999. epsi=1. cdie v14f=8. e14f=2. action ~ ~ To run this input type (once you have a path in your .cshrc leading to moil.src.001/exe, and assuming you called the file "input" again): energy output this will activate the part of MOIL that deals with the energy calculation. The output looks like this (note the particle pairs are counted for the non-bonded interactions): -------------------------------------------------------------------------- pote> file rcon name=(BENZ.WCON) unit=10 read pote> file rcrd name=(BENZ.CRD) unit=11 read pote> file wene name=(BENZ.ENE) unit=13 wovr pote> rmax=9999. epsi=1. cdie pote> action getcrd> * benzene getcrd> * Number of particle pairs 0 ----------------------------------------------------------------------- The results of the energy calculation are in the file BENZ.ENE : --------------------------------------------------------------------------- Parameter for energy calculation Constant dielectric will be used. Cutoff= 9999.00000 Dielectric constant= 1.00000 1-4 scaling factors van der Waals= 1.00000 Electrostatic= 1.00000 ENERGIES: E total = 1.012 E bond = .991 E angl = .006 E tors = .014 E impr = .000 E vdw = .000 E elec = .000 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = 10.197 --------------------------------------------------------------------------- In order to minimize this energy, you activate the MOIL minimizer by typing: mini output This time the input looks like this: ~ ~ input for energy minimization ~ file rcon name=(BENZ.WCON) unit=10 read file rcrd name=(BENZ.CRD) unit=11 read file wcrd name=(BENZMIN.CRD) unit=12 wovr file wmin name=(BENZ.MIN) unit=13 wovr rmax=9999. epsi=50. cdie tolf=0.01 mistep=200 list=50 action ~ ~ The connectivity and coordinate files make up your molecule; the minimization history is written in *.MIN and the energy minimized coordinates are written in *MIN.CRD Here is the output you get: -------------------------------------------------------------------------- mini> file rcon name=(BENZ.WCON) unit=10 read mini> file rcrd name=(BENZ.CRD) unit=11 read mini> file wcrd name=(BENZMIN.CRD) unit=12 wovr mini> file wmin name=(BENZ.MIN) unit=13 wovr mini> rmax=9999. epsi=50. cdie mini> tolf=0.01 mistep=200 list=50 mini> action Minimization parameters Maximum number of minimization steps 200 Update non-bond list each 50 steps Connectivity file to be read from unit 10 Initial coordinates to be read from unit 11 Final coordinates will be written to unit 12 The cutoff distance (for elec. and vdw) is 9999.00000 The second cutoff (for neighbor buffer) is ********** the allowed gradient error is .01000 getcrd> * benzene getcrd> * Number of particle pairs 0 *** Minimization completed ---------------------------------------------------------------------- And here are the results of the minimization: a) file BENZ.MIN ---------------------------------------------------------------------- At step 0 Parameter for energy calculation Constant dielectric will be used. Cutoff= 9999.00000 Dielectric constant= 50.00000 1-4 scaling factors van der Waals= 1.00000 Electrostatic= 1.00000 After 21 energy calls Current gradient norm .00193 requested .01000 Minimization converged ENERGIES: E total = .000 E bond = .000 E angl = .000 E tors = .000 E impr = .000 E vdw = .000 E elec = .000 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = .002 ---------------------------------------------------------------------- b) file BENZMIN.CRD (the energy minimized coordinates) --------------------------------------------------------------------- * title for CHARMM coordinates * 6 1 1 BENZ C1 10.10735 46.18345 12.15675 BENZ FREE 13.01900 2 1 BENZ C2 11.13765 46.27098 11.21292 BENZ FREE 13.01900 3 1 BENZ C3 10.97774 47.06940 10.07408 BENZ FREE 13.01900 4 1 BENZ C4 9.78753 47.78030 9.87909 BENZ FREE 13.01900 5 1 BENZ C5 8.75723 47.69276 10.82291 BENZ FREE 13.01900 6 1 BENZ C6 8.91713 46.89434 11.96175 BENZ FREE 13.01900 ----------------------------------------------------------------------- This is all we shall do with benzene. The next example will move us closer to working with proteins, which is what MOIL is really about. ======================================================================= ======================================================================= Let's start with a dipeptide. This time we will generate the molecule using the database that came with MOIL. ALL.MONO lists the monomer amino-acids and various kinds of termini. ALL_PROP lists the properties of atoms. As before we start with setting up the connectivity for the molecule. Here is the input: ~ ~ input for connectivity ~ file prop name=(ALL_PROP) unit=10 read file mono name=(ALL.MONO) unit=11 read file poly name=(VALD.POLY) unit=12 read file wcon name=(VALD.WCON) unit=13 wovr action *EOD ~ ~ The files ALL_PROP and ALL.MONO you already have. Here is the polymerization file VALD.POLY where you tell MOIL how many monomers you are joining and in what sequence. ---------------------------------------------------------------------- ~ MOLC=(VALD) #mon=3 NTR0 VAL CTR0 *EOD ~ ---------------------------------------------------------------------- Each end of the amino-acid has a neutral terminus attached to it (the termini are also defined as monomers). All three joined together will have 2 peptide bonds, hence the valine dipeptide, if you weren't sure. You can find the names of other termini in ALL.MONO whenever you need them. You run the above input just like before (activate the connectivity subprogram of MOIL): con output This long output simply echoes the contents of the data files: ----------------------------------------------------------------------- CONN> file prop name=(ALL_PROP) unit=10 read CONN> file mono name=(ALL.MONO) unit=11 read CONN> file poly name=(VALD.POLY) unit=12 read CONN> file wcon name=(VALD.WCON) unit=13 wovr CONN> action propert> PRTC . . . propert> DONE *** DONE WITH IMPROPER TORSIONS *** *** Number of parameters read: Unique particles, bonds, angles, torsions, improper torsions 100 155 351 344 163 monomers> MONO LIST monomers> ~ put in ALL.MONO file monomers> ~ monomers> ~uncharged N-terminus monomers> ~ monomers> MONO=(NTR0) #prt=6 chrg=0. monomers> ~ monomers> ~ X monomers> ~ : : monomers> ~ CH3-CO - NH.CAH .. C..O monomers> ~ || : monomers> ~ OC HN monomers> ~ monomers> UNIQ=(N) PRTC=(NH) NEXT monomers> UNIQ=(CA) PRTC=(CAH) NEXT monomers> UNIQ=(H) PRTC=(HN) NEXT monomers> UNIQ=(ME) PRTC=(CH3) monomers> UNIQ=(C) PRTC=(CO) monomers> UNIQ=(O) PRTC=(OC) monomers> DONE monomers> BOND monomers> C-N* C-O C-ME monomers> DONE monomers> ~ monomers> *EOD LEGO> MOLC=(VALD) #mon=3 LEGO> NTR0 VAL CTR0 *** MOLECULE VALD ASSEMBLED *** LEGO> *EOD ----------------------------------------------------------------------- Here is the connectivity file that results from all this reading: ---------------------------------------------------------------------- ~ CONNECTIVITY FILE FOR MOLECULES: ~ totmon npt nb nangl ntors nimp totex totspe lestyp NBULK 3 14 13 18 13 6 31 20 0 1 VALD ~ Pointers to last particle of BULK 14 ~ Monomer names NTR0 VAL CTR0 ~ Pointers to last particle of monomer 3 11 14 ~ Properties of particles list : ~pt mono ptid lesid ptnm ptms ptchg epsgm6 epsgm12 ptwei 1 1 13 0 ME 15.00 .00000 .47821E+02 .28586E+04 .10000E+01 2 1 9 0 C 12.00 .50000 .34176E+02 .18022E+04 .10000E+01 3 1 10 0 O 16.00 -.50000 .23769E+02 .61644E+03 .10000E+01 4 2 7 0 N 14.00 -.57000 .28308E+02 .97175E+03 .10000E+01 5 2 8 0 H 1.00 .37000 .00000E+00 .00000E+00 .10000E+01 6 2 11 0 CA 13.00 .20000 .31040E+02 .17032E+04 .10000E+01 7 2 49 0 CB 13.00 .00000 .32282E+02 .18422E+04 .10000E+01 8 2 26 0 CG1 15.00 .00000 .47821E+02 .28586E+04 .10000E+01 9 2 13 0 CG2 15.00 .00000 .47821E+02 .28586E+04 .10000E+01 10 2 9 0 C 12.00 .50000 .34176E+02 .18022E+04 .10000E+01 11 2 10 0 O 16.00 -.50000 .23769E+02 .61644E+03 .10000E+01 12 3 7 0 N 14.00 -.57000 .28308E+02 .97175E+03 .10000E+01 13 3 8 0 H 1.00 .37000 .00000E+00 .00000E+00 .10000E+01 14 3 12 0 ME 15.00 .20000 .45249E+02 .24829E+04 .10000E+01 ~ Bonds list: ~ ib1 ib2 kbond req 1 2 317.0000 1.5220 2 3 317.0000 1.2290 2 4 317.0000 1.3350 4 6 260.0000 1.4490 4 5 337.0000 1.0100 6 7 260.0000 1.5260 6 10 317.0000 1.5220 7 9 260.0000 1.5260 7 8 260.0000 1.5260 10 12 490.0000 1.3350 10 11 570.0000 1.2290 12 14 337.0000 1.4490 12 13 434.0000 1.0100 ~ Angles list: ~ iangl1 iangl2 iangl3 kangl angleq 1 2 3 80.00000 120.40000 1 2 4 50.00000 121.90000 3 2 4 80.00000 122.90000 2 4 5 35.00000 119.80000 2 4 6 50.00000 121.90000 5 4 6 38.00000 118.40000 4 6 7 80.00000 109.70000 4 6 10 63.00000 110.10000 7 6 10 63.00000 111.10000 6 7 8 63.00000 111.50000 6 7 9 63.00000 111.50000 8 7 9 63.00000 111.50000 6 10 11 80.00000 120.40000 6 10 12 70.00000 116.60000 11 10 12 80.00000 122.90000 10 12 13 35.00000 119.80000 10 12 14 50.00000 121.90000 13 12 14 38.00000 118.40000 ~ Torsions list: ~ itor1 itor2 itor3 itor4 period ktors1 ktors2 ktors3 phase 1 2 4 5 2 .0000 2.5000 .0000 -1.000 1 2 4 6 2 .0000 2.5000 .0000 -1.000 3 2 4 5 2 .0000 2.5000 .0000 -1.000 3 2 4 6 2 .0000 2.5000 .0000 -1.000 4 6 7 8 3 .0000 .0000 .5000 -1.000 4 6 7 9 3 .0000 .0000 .5000 1.000 4 6 10 11 3 .0000 .0000 .1000 -1.000 8 7 6 10 3 .0000 .0000 .5000 1.000 9 7 6 10 3 .0000 .0000 .5000 1.000 6 10 12 13 2 .0000 2.5000 .0000 -1.000 6 10 12 14 2 .0000 2.5000 .0000 -1.000 11 10 12 13 2 .0000 2.5000 .0000 -1.000 11 10 12 14 2 .0000 2.5000 .0000 -1.000 ~ Improper torsion properties: ~iimp1 iimp2 iimp3 iimp4 kimp impeq 2 1 3 4 .10000000E+03 .00000000E+00 4 6 2 5 .45000000E+02 .00000000E+00 6 4 10 7 .55000000E+02 .35260000E+02 7 6 9 8 .55000000E+02 .35260000E+02 10 6 11 12 .10000000E+03 .00000000E+00 12 14 10 13 .45000000E+02 .00000000E+00 ~ Exclusion list 1-2 1-3, set as followed: ~ atom number, number of exclusions and list 1 3 2 3 4 2 4 3 4 5 6 3 1 4 4 4 6 5 7 10 5 1 6 6 6 7 10 8 9 11 12 7 3 9 8 10 8 1 9 10 4 12 11 13 14 11 1 12 12 2 14 13 13 1 14 ~ Special list 1-4 set as followed: ~ atom number, number of exclusions and list 1 2 5 6 2 2 7 10 3 2 5 6 4 4 8 9 11 12 5 2 7 10 6 2 13 14 7 2 11 12 8 1 10 9 1 10 11 2 13 14 ----------------------------------------------------------------------- The connectivity file must be supplemented by the atomic coordinates' file of course. What's the use of having the connections if you cannot make some-body utilize them? Here are the coordinates I have: VH.CRD ----------------------------------------------------------------------- * initial structure, valine dipeptide in helix conformation * 14 1 1 NTR0 ME -0.02717 3.41564 0.00488 VALD 1 15.03500 2 1 NTR0 C 0.04827 1.91708 -0.24090 VALD 1 12.01100 3 1 NTR0 O 0.67142 1.46396 -1.20000 VALD 1 15.99940 4 2 VAL N -0.58751 1.13019 0.63714 VALD 1 14.00670 5 2 VAL H -1.05774 1.57408 1.37031 VALD 1 1.00800 6 2 VAL CA -0.63665 -0.33703 0.58146 VALD 1 13.01900 7 2 VAL CB -1.49790 -0.83846 -0.61850 VALD 1 15.03500 8 2 VAL CG1 -2.91731 -0.28641 -0.56404 VALD 1 0.00000 9 2 VAL CG2 -1.57439 -2.35741 -0.74809 VALD 1 0.00000 10 2 VAL C 0.70303 -1.08091 0.68144 VALD 1 12.01100 11 2 VAL O 0.92822 -1.80836 1.64793 VALD 1 15.99940 12 3 CTR0 N 1.61106 -0.91246 -0.28244 VALD 1 14.00670 13 3 CTR0 H 1.39819 -0.28295 -1.00480 VALD 1 1.00800 14 3 CTR0 ME 2.93848 -1.59696 -0.26438 VALD 1 15.03500 ----------------------------------------------------------------------- This is a valine dipeptide in the right-handed helix conformation. Here are the coordinates of the same dipeptide, but in a beta-sheet conformation. We shall need both of them for the next example. VB.CRD ----------------------------------------------------------------------- * final structure, valine dipeptide in beta-sheet conformation * 14 1 1 NTR0 ME 0.03442 3.38506 -0.46052 VALD 1 15.03500 2 1 NTR0 C 0.18913 2.00194 0.15123 VALD 1 12.01100 3 1 NTR0 O 1.04729 1.77335 1.00482 VALD 1 15.99940 4 2 VAL N -0.60845 1.03091 -0.30180 VALD 1 14.00670 5 2 VAL H -1.24762 1.24171 -1.01118 VALD 1 1.00800 6 2 VAL CA -0.53889 -0.33325 0.22054 VALD 1 13.01900 7 2 VAL CB -1.89084 -1.06501 -0.00613 VALD 1 15.03500 8 2 VAL CG1 -3.07008 -0.27380 0.54780 VALD 1 0.00000 9 2 VAL CG2 -1.92494 -2.46055 0.61079 VALD 1 0.00000 10 2 VAL C 0.66176 -1.10417 -0.34834 VALD 1 12.01100 11 2 VAL O 0.56886 -1.89159 -1.29147 VALD 1 15.99940 12 3 CTR0 N 1.83107 -0.79657 0.21331 VALD 1 14.00670 13 3 CTR0 H 1.83535 -0.07016 0.87597 VALD 1 1.00800 14 3 CTR0 ME 3.11294 -1.43787 -0.20500 VALD 1 15.03500 ----------------------------------------------------------------------- This time we will make MOIL work a little harder - we will calculate the reaction path for the conformational change between the right-handed helix and the beta-sheet using the two valine dipeptide conformations. For this we need to place both conformations in their respective minima, and then determine the path they take on the potential energy map as they move from one conformation to the other. First, let's see if the structures are reasonable, i.e., let's check their energies. The input goes like this: ~ ~ input - energy test for the initial structure (helix): VH.CRD ~ file rcon name=(VALD.WCON) unit=10 read file rcrd name=(VH.CRD) unit=11 read file wene name=(VH.ENE) unit=13 wovr rmax=9999. epsi=1. cdie v14f=8. e14f=2. action ~ Run it as usual: energy < input > output And the output file VH.ENE gives: -------------------------------------------------------------------------- Parameter for energy calculation Constant dielectric will be used. Cutoff= 9999.00000 Dielectric constant= 1.00000 1-4 scaling factors van der Waals= 8.00000 Electrostatic= 2.00000 ENERGIES: E total = 35.767 E bond = 1.848 E angl = 2.303 E tors = 1.230 E impr = 78.916 E vdw = 3.996 E elec = -52.527 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = 31.318 -------------------------------------------------------------------------- Similarly for the final structure on the path: VB.CRD ~ ~ input - energy test for beta sheet conformation ~ file rcon name=(VALD.WCON) unit=10 read file rcrd name=(VH.CRD) unit=11 read file wene name=(VH.ENE) unit=13 wovr rmax=9999. epsi=1. cdie v14f=8. e14f=2. action ~ and the output is: --------------------------------------------------------------------------- Parameter for energy calculation Constant dielectric will be used. Cutoff= 9999.00000 Dielectric constant= 1.00000 1-4 scaling factors van der Waals= 8.00000 Electrostatic= 2.00000 ENERGIES: E total = 34.489 E bond = 1.638 E angl = 1.191 E tors = 1.215 E impr = 79.806 E vdw = 2.988 E elec = -52.349 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = 31.490 ------------------------------------------------------------------------- Now you see that these two structures are not quite at the minima on the potential energy surface that want them to be sitting at. Here is the input we shall use to minimize their energy: ~ ~ minimization of helix conformation ~ ~debug file rcon name=(VALD.WCON) unit=10 read file rcrd name=(VH.CRD) unit=11 read file wcrd name=(VHMIN.CRD) unit=12 wovr file wmin name=(VH.MIN) unit=13 wovr rmax=9999. epsi=1. cdie v14f=8. e14f=2. tolf=0.0001 mistep=1100 list=50 action ~ ~ (again, the description of the variables can be found in file mini.prog) To run it type (substituting your file names, of course): mini < input >output The output looks like this: --------------------------------------------------------------------------- mini> ~ minimization of helix conformation mini> ~ mini> ~debug mini> file rcon name=(VALD.WCON) unit=10 read mini> file rcrd name=(VH.CRD) unit=11 read mini> file wcrd name=(VHMIN.CRD) unit=12 wovr mini> file wmin name=(VH.MIN) unit=13 wovr mini> rmax=9999. epsi=1. cdie v14f=8. e14f=2. mini> tolf=0.0001 mistep=1100 list=50 mini> action Minimization parameters Maximum number of minimization steps 1100 Update non-bond list each 50 steps Connectivity file to be read from unit 10 Initial coordinates to be read from unit 11 Final coordinates will be written to unit 12 The cutoff distance (for elec. and vdw) is 9999.00000 The second cutoff (for neighbor buffer) is ********** the allowed gradient error is .00010 getcrd> * initial structure, valine dipeptide in helix conformation getcrd> * Number of particle pairs 40 Number of particle pairs 40 Number of particle pairs 40 Number of particle pairs 40 Number of particle pairs 40 Number of particle pairs 40 Number of particle pairs 40 *** Minimization completed ----------------------------------------------------------------------- And here is the file we wrote the minimization history into (VH.MIN): ------------------------------------------------------------------------ At step 0 Parameter for energy calculation Constant dielectric will be used. Cutoff= 9999.00000 Dielectric constant= 1.00000 1-4 scaling factors van der Waals= 8.00000 Electrostatic= 2.00000 After 50 energy calls Current gradient norm 3.12341 requested .00010 ENERGIES: E total = -38.339 E bond = 1.160 E angl = 4.320 E tors = 2.442 E impr = 2.156 E vdw = 6.382 E elec = -54.798 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = 3.123 After 100 energy calls Current gradient norm 1.17786 requested .00010 ENERGIES: E total = -42.358 E bond = .586 E angl = 4.127 E tors = 1.234 E impr = 1.437 E vdw = 4.229 E elec = -53.971 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = 1.178 After 150 energy calls Current gradient norm .12264 requested .00010 ENERGIES: E total = -42.546 E bond = .563 E angl = 4.234 E tors = 1.204 E impr = 1.361 E vdw = 4.254 E elec = -54.163 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = .123 After 200 energy calls Current gradient norm .01096 requested .00010 ENERGIES: E total = -42.552 E bond = .556 E angl = 4.293 E tors = 1.203 E impr = 1.352 E vdw = 4.148 E elec = -54.105 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = .011 After 250 energy calls Current gradient norm .00134 requested .00010 ENERGIES: E total = -42.552 E bond = .557 E angl = 4.296 E tors = 1.203 E impr = 1.351 E vdw = 4.149 E elec = -54.108 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = .001 After 300 energy calls Current gradient norm .00026 requested .00010 ENERGIES: E total = -42.552 E bond = .557 E angl = 4.297 E tors = 1.203 E impr = 1.350 E vdw = 4.148 E elec = -54.107 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = .000 After 350 energy calls Current gradient norm .00004 requested .00010 Minimization converged ENERGIES: E total = -42.552 E bond = .557 E angl = 4.297 E tors = 1.203 E impr = 1.350 E vdw = 4.148 E elec = -54.107 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = .000 -------------------------------------------------------------------------- Now for the beta-sheet structure. This time we will write the energies less often (list = "every 100 steps"). ~ ~ input - minimization of the beta-sheet conformation ~ ~debug file rcon name=(VALD.WCON) unit=10 read file rcrd name=(VB.CRD) unit=11 read file wcrd name=(VBMIN.CRD) unit=12 wovr file wmin name=(VB.MIN) unit=13 wovr rmax=9999. epsi=1. cdie v14f=8. e14f=2. tolf=0.0001 mistep=1100 list=100 action ~ and output: --------------------------------------------------------------------------- mini> ~ input - minimization of beta-sheet conformation mini> ~ mini> ~debug mini> file rcon name=(VALD.WCON) unit=10 read mini> file rcrd name=(VH.CRD) unit=11 read mini> file wcrd name=(VHMIN.CRD) unit=12 wovr mini> file wmin name=(VH.MIN) unit=13 wovr mini> rmax=9999. epsi=1. cdie v14f=8. e14f=2. mini> tolf=0.0001 mistep=1100 list=50 mini> action Minimization parameters Maximum number of minimization steps 1100 Update non-bond list each 50 steps Connectivity file to be read from unit 10 Initial coordinates to be read from unit 11 Final coordinates will be written to unit 12 The cutoff distance (for elec. and vdw) is 9999.00000 The second cutoff (for neighbor buffer) is ********** the allowed gradient error is .00010 getcrd> * initial structure, valine dipeptide in helix conformation getcrd> * Number of particle pairs 40 Number of particle pairs 40 Number of particle pairs 40 Number of particle pairs 40 Number of particle pairs 40 Number of particle pairs 40 Number of particle pairs 40 *** Minimization completed -------------------------------------------------------------------------- And the minimization history (VB.CRD): --------------------------------------------------------------------------- At step 0 Parameter for energy calculation Constant dielectric will be used. Cutoff= 9999.00000 Dielectric constant= 1.00000 1-4 scaling factors van der Waals= 8.00000 Electrostatic= 2.00000 After 100 energy calls Current gradient norm .52888 requested .00010 ENERGIES: E total = -45.183 E bond = .251 E angl = 1.317 E tors = 1.267 E impr = .945 E vdw = 1.645 E elec = -50.607 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = .529 After 200 energy calls Current gradient norm .19795 requested .00010 ENERGIES: E total = -45.846 E bond = .325 E angl = 1.619 E tors = 1.176 E impr = .807 E vdw = 3.055 E elec = -52.829 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = .198 After 300 energy calls Current gradient norm .00300 requested .00010 ENERGIES: E total = -45.851 E bond = .302 E angl = 1.621 E tors = 1.180 E impr = .815 E vdw = 2.950 E elec = -52.719 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = .003 After 400 energy calls Current gradient norm .00064 requested .00010 ENERGIES: E total = -45.851 E bond = .302 E angl = 1.620 E tors = 1.181 E impr = .814 E vdw = 2.952 E elec = -52.721 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = .001 After 449 energy calls Current gradient norm .00001 requested .00010 Minimization converged ENERGIES: E total = -45.851 E bond = .302 E angl = 1.620 E tors = 1.181 E impr = .814 E vdw = 2.952 E elec = -52.721 E cnst = .000 E evsym= .000 E elsym= .000 Norm Force = .000 ------------------------------------------------------------------------ We finally have the initial (reactant) and final (product) structures where we want them, so the time has come for the path calculation. ========================================================================= CHMIN (minimum energy reaction path) ========================================================================= We will find the minimum energy path for the transformation between the "reactants" and the "products" using the SPW (Self Penalty Walk) chain of Czerminski and Elber (Int.J.Quant.Chem. 24, 167, 1990). First the program will interpolate between the initial and the final structures by constructing a set of intermediary structures. (you can also look at the structures as monomers in a linear polymer) These "chain" structures are "intraconnected" by a equidistance constraints and repulsions between neighbors. Then the energy of the chain will be minimized and the resulting structures lying along the minimum energy path will be written into a binary file. You can find the description of the path format, and of the parameters involved, in the file chain.prog (it describes a different variant of the same calculation). Now let's do the calculation. The input is: ~ ~ input f