183da9db09c0cc7fa61d69c1b08b878b62575ae8
[unres.git] / parmread.F
1       subroutine parmread
2 C
3 C Read the parameters of the probability distributions of the virtual-bond
4 C valence angles and the side chains and energy parameters.
5 C
6 C Important! Energy-term weights ARE NOT read here; they are read from the
7 C main input file instead, because NO defaults have yet been set for these
8 C parameters.
9 C
10       implicit real*8 (a-h,o-z)
11       include 'DIMENSIONS'
12 #ifdef MPI
13       include "mpif.h"
14       integer IERROR
15 #endif
16       include 'COMMON.IOUNITS'
17       include 'COMMON.CHAIN'
18       include 'COMMON.INTERACT'
19       include 'COMMON.GEO'
20       include 'COMMON.LOCAL'
21       include 'COMMON.TORSION'
22       include 'COMMON.SCCOR'
23       include 'COMMON.SCROT'
24       include 'COMMON.FFIELD'
25       include 'COMMON.NAMES'
26       include 'COMMON.SBRIDGE'
27       include 'COMMON.MD'
28       include 'COMMON.SETUP'
29       character*1 t1,t2,t3
30       character*1 onelett(4) /"G","A","P","D"/
31       logical lprint,LaTeX
32       dimension blower(3,3,maxlob)
33       dimension b(13)
34       character*3 lancuch,ucase
35 C
36 C For printing parameters after they are read set the following in the UNRES
37 C C-shell script:
38 C
39 C setenv PRINT_PARM YES
40 C
41 C To print parameters in LaTeX format rather than as ASCII tables:
42 C
43 C setenv LATEX YES
44 C
45       call getenv_loc("PRINT_PARM",lancuch)
46       lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
47       call getenv_loc("LATEX",lancuch)
48       LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
49 C
50       dwa16=2.0d0**(1.0d0/6.0d0)
51       itypro=20
52 C Assign virtual-bond length
53       vbl=3.8D0
54       vblinv=1.0D0/vbl
55       vblinv2=vblinv*vblinv
56 c
57 c Read the virtual-bond parameters, masses, and moments of inertia
58 c and Stokes' radii of the peptide group and side chains
59 c
60 #ifdef CRYST_BOND
61       read (ibond,*) vbldp0,akp,mp,ip,pstok
62       do i=1,ntyp
63         nbondterm(i)=1
64         read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
65         dsc(i) = vbldsc0(1,i)
66         if (i.eq.10) then
67           dsc_inv(i)=0.0D0
68         else
69           dsc_inv(i)=1.0D0/dsc(i)
70         endif
71       enddo
72 #else
73       read (ibond,*) junk,vbldp0,akp,rjunk,mp,ip,pstok
74       do i=1,ntyp
75         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
76      &   j=1,nbondterm(i)),msc(i),isc(i),restok(i)
77         dsc(i) = vbldsc0(1,i)
78         if (i.eq.10) then
79           dsc_inv(i)=0.0D0
80         else
81           dsc_inv(i)=1.0D0/dsc(i)
82         endif
83       enddo
84 #endif
85       if (lprint) then
86         write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
87         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',
88      &   'inertia','Pstok'
89         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp,ip,pstok
90         do i=1,ntyp
91           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
92      &      vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i),isc(i),restok(i)
93           do j=2,nbondterm(i)
94             write (iout,'(13x,3f10.5)')
95      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
96           enddo
97         enddo
98       endif
99 #ifdef CRYST_THETA
100 C
101 C Read the parameters of the probability distribution/energy expression 
102 C of the virtual-bond valence angles theta
103 C
104       do i=1,ntyp
105         read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
106      &    (bthet(j,i),j=1,2)
107         read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
108         read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
109         read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
110         sigc0(i)=sigc0(i)**2
111       enddo
112       close (ithep)
113       if (lprint) then
114       if (.not.LaTeX) then
115         write (iout,'(a)') 
116      &    'Parameters of the virtual-bond valence angles:'
117         write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
118      & '    ATHETA0   ','         A1   ','        A2    ',
119      & '        B1    ','         B2   '        
120         do i=1,ntyp
121           write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
122      &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
123         enddo
124         write (iout,'(/a/9x,5a/79(1h-))') 
125      & 'Parameters of the expression for sigma(theta_c):',
126      & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
127      & '     ALPH3    ','    SIGMA0C   '        
128         do i=1,ntyp
129           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
130      &      (polthet(j,i),j=0,3),sigc0(i) 
131         enddo
132         write (iout,'(/a/9x,5a/79(1h-))') 
133      & 'Parameters of the second gaussian:',
134      & '    THETA0    ','     SIGMA0   ','        G1    ',
135      & '        G2    ','         G3   '        
136         do i=1,ntyp
137           write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
138      &       sig0(i),(gthet(j,i),j=1,3)
139         enddo
140        else
141         write (iout,'(a)') 
142      &    'Parameters of the virtual-bond valence angles:'
143         write (iout,'(/a/9x,5a/79(1h-))') 
144      & 'Coefficients of expansion',
145      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
146      & '   b1*10^1    ','    b2*10^1   '        
147         do i=1,ntyp
148           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
149      &        a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
150         enddo
151         write (iout,'(/a/9x,5a/79(1h-))') 
152      & 'Parameters of the expression for sigma(theta_c):',
153      & ' alpha0       ','  alph1       ',' alph2        ',
154      & ' alhp3        ','   sigma0c    '        
155         do i=1,ntyp
156           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
157      &      (polthet(j,i),j=0,3),sigc0(i) 
158         enddo
159         write (iout,'(/a/9x,5a/79(1h-))') 
160      & 'Parameters of the second gaussian:',
161      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
162      & '        G2    ','   G3*10^1    '        
163         do i=1,ntyp
164           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
165      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
166         enddo
167       endif
168       endif
169 #else 
170
171 C Read the parameters of Utheta determined from ab initio surfaces
172 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
173 C
174       read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,
175      &  ntheterm3,nsingle,ndouble
176       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
177       read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
178       do i=1,maxthetyp
179         do j=1,maxthetyp
180           do k=1,maxthetyp
181             aa0thet(i,j,k)=0.0d0
182             do l=1,ntheterm
183               aathet(l,i,j,k)=0.0d0
184             enddo
185             do l=1,ntheterm2
186               do m=1,nsingle
187                 bbthet(m,l,i,j,k)=0.0d0
188                 ccthet(m,l,i,j,k)=0.0d0
189                 ddthet(m,l,i,j,k)=0.0d0
190                 eethet(m,l,i,j,k)=0.0d0
191               enddo
192             enddo
193             do l=1,ntheterm3
194               do m=1,ndouble
195                 do mm=1,ndouble
196                  ffthet(mm,m,l,i,j,k)=0.0d0
197                  ggthet(mm,m,l,i,j,k)=0.0d0
198                 enddo
199               enddo
200             enddo
201           enddo
202         enddo
203       enddo 
204       do i=1,nthetyp
205         do j=1,nthetyp
206           do k=1,nthetyp
207             read (ithep,'(3a)',end=111,err=111) res1,res2,res3
208             read (ithep,*,end=111,err=111) aa0thet(i,j,k)
209             read (ithep,*,end=111,err=111)(aathet(l,i,j,k),l=1,ntheterm)
210             read (ithep,*,end=111,err=111)
211      &       ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
212      &        (ccthet(lll,ll,i,j,k),lll=1,nsingle),
213      &        (ddthet(lll,ll,i,j,k),lll=1,nsingle),
214      &        (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
215             read (ithep,*,end=111,err=111)
216      &      (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
217      &         ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
218      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
219           enddo
220         enddo
221       enddo
222 C
223 C For dummy ends assign glycine-type coefficients of theta-only terms; the
224 C coefficients of theta-and-gamma-dependent terms are zero.
225 C
226       do i=1,nthetyp
227         do j=1,nthetyp
228           do l=1,ntheterm
229             aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
230             aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
231           enddo
232           aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
233           aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
234         enddo
235         do l=1,ntheterm
236           aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
237         enddo
238         aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
239       enddo
240 C
241 C Control printout of the coefficients of virtual-bond-angle potentials
242 C
243       if (lprint) then
244         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
245         do i=1,nthetyp+1
246           do j=1,nthetyp+1
247             do k=1,nthetyp+1
248               write (iout,'(//4a)') 
249      &         'Type ',onelett(i),onelett(j),onelett(k) 
250               write (iout,'(//a,10x,a)') " l","a[l]"
251               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
252               write (iout,'(i2,1pe15.5)')
253      &           (l,aathet(l,i,j,k),l=1,ntheterm)
254             do l=1,ntheterm2
255               write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') 
256      &          "b",l,"c",l,"d",l,"e",l
257               do m=1,nsingle
258                 write (iout,'(i2,4(1pe15.5))') m,
259      &          bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
260      &          ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
261               enddo
262             enddo
263             do l=1,ntheterm3
264               write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))')
265      &          "f+",l,"f-",l,"g+",l,"g-",l
266               do m=2,ndouble
267                 do n=1,m-1
268                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
269      &              ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
270      &              ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
271                 enddo
272               enddo
273             enddo
274           enddo
275         enddo
276       enddo
277       call flush(iout)
278       endif
279       write (2,*) "Start reading THETA_PDB"
280       do i=1,ntyp
281         read (ithep_pdb,*,err=111,end=111) a0thet(i),(athet(j,i),j=1,2),
282      &    (bthet(j,i),j=1,2)
283         read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
284         read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
285         read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
286         sigc0(i)=sigc0(i)**2
287       enddo
288       write (2,*) "End reading THETA_PDB"
289       close (ithep_pdb)
290 #endif
291       close(ithep)
292 #ifdef CRYST_SC
293 C
294 C Read the parameters of the probability distribution/energy expression
295 C of the side chains.
296 C
297       do i=1,ntyp
298         read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
299         if (i.eq.10) then
300           dsc_inv(i)=0.0D0
301         else
302           dsc_inv(i)=1.0D0/dsc(i)
303         endif
304         if (i.ne.10) then
305         do j=1,nlob(i)
306           do k=1,3
307             do l=1,3
308               blower(l,k,j)=0.0D0
309             enddo
310           enddo
311         enddo  
312         bsc(1,i)=0.0D0
313         read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),
314      &    ((blower(k,l,1),l=1,k),k=1,3)
315         do j=2,nlob(i)
316           read (irotam,*,end=112,err=112) bsc(j,i)
317           read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),
318      &                                 ((blower(k,l,j),l=1,k),k=1,3)
319         enddo
320         do j=1,nlob(i)
321           do k=1,3
322             do l=1,k
323               akl=0.0D0
324               do m=1,3
325                 akl=akl+blower(k,m,j)*blower(l,m,j)
326               enddo
327               gaussc(k,l,j,i)=akl
328               gaussc(l,k,j,i)=akl
329             enddo
330           enddo 
331         enddo
332         endif
333       enddo
334       close (irotam)
335       if (lprint) then
336         write (iout,'(/a)') 'Parameters of side-chain local geometry'
337         do i=1,ntyp
338           nlobi=nlob(i)
339           if (nlobi.gt.0) then
340             if (LaTeX) then
341               write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
342      &         ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
343                write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
344      &                             'log h',(bsc(j,i),j=1,nlobi)
345                write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
346      &        'x',((censc(k,j,i),k=1,3),j=1,nlobi)
347               do k=1,3
348                 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
349      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
350               enddo
351             else
352               write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
353               write (iout,'(a,f10.4,4(16x,f10.4))')
354      &                             'Center  ',(bsc(j,i),j=1,nlobi)
355               write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),
356      &           j=1,nlobi)
357               write (iout,'(a)')
358             endif
359           endif
360         enddo
361       endif
362 #else
363
364 C Read scrot parameters for potentials determined from all-atom AM1 calculations
365 C added by Urszula Kozlowska 07/11/2007
366 C
367       do i=1,ntyp
368         read (irotam,*,end=112,err=112) 
369        if (i.eq.10) then 
370          read (irotam,*,end=112,err=112) 
371        else
372          do j=1,65
373            read(irotam,*,end=112,err=112) sc_parmin(j,i)
374          enddo  
375        endif
376       enddo
377 C
378 C Read the parameters of the probability distribution/energy expression
379 C of the side chains.
380 C
381       do i=1,ntyp
382         read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
383         if (i.eq.10) then
384           dsc_inv(i)=0.0D0
385         else
386           dsc_inv(i)=1.0D0/dsc(i)
387         endif
388         if (i.ne.10) then
389         do j=1,nlob(i)
390           do k=1,3
391             do l=1,3
392               blower(l,k,j)=0.0D0
393             enddo
394           enddo
395         enddo  
396         bsc(1,i)=0.0D0
397         read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),
398      &    ((blower(k,l,1),l=1,k),k=1,3)
399         do j=2,nlob(i)
400           read (irotam_pdb,*,end=112,err=112) bsc(j,i)
401           read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),
402      &                                 ((blower(k,l,j),l=1,k),k=1,3)
403         enddo
404         do j=1,nlob(i)
405           do k=1,3
406             do l=1,k
407               akl=0.0D0
408               do m=1,3
409                 akl=akl+blower(k,m,j)*blower(l,m,j)
410               enddo
411               gaussc(k,l,j,i)=akl
412               gaussc(l,k,j,i)=akl
413             enddo
414           enddo 
415         enddo
416         endif
417       enddo
418       close (irotam_pdb)
419 #endif
420       close(irotam)
421
422 #ifdef CRYST_TOR
423 C
424 C Read torsional parameters in old format
425 C
426       read (itorp,*,end=113,err=113) ntortyp,nterm_old
427       if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
428       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
429       do i=1,ntortyp
430         do j=1,ntortyp
431           read (itorp,'(a)')
432           do k=1,nterm_old
433             read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i) 
434           enddo
435         enddo
436       enddo
437       close (itorp)
438       if (lprint) then
439         write (iout,'(/a/)') 'Torsional constants:'
440         do i=1,ntortyp
441           do j=1,ntortyp
442             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
443             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
444           enddo
445         enddo
446       endif
447 #else
448 C
449 C Read torsional parameters
450 C
451       read (itorp,*,end=113,err=113) ntortyp
452       read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
453       itortyp(ntyp1)=ntortyp+1
454 c      write (iout,*) 'ntortyp',ntortyp
455       do i=1,ntortyp
456         do j=1,ntortyp
457           read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j)
458           v0ij=0.0d0
459           si=-1.0d0
460           do k=1,nterm(i,j)
461             read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j) 
462             v0ij=v0ij+si*v1(k,i,j)
463             si=-si
464           enddo
465           do k=1,nlor(i,j)
466             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
467      &        vlor2(k,i,j),vlor3(k,i,j) 
468             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
469           enddo
470           v0(i,j)=v0ij
471         enddo
472       enddo
473       close (itorp)
474       if (lprint) then
475         write (iout,'(/a/)') 'Torsional constants:'
476         do i=1,ntortyp
477           do j=1,ntortyp
478             write (iout,*) 'ityp',i,' jtyp',j
479             write (iout,*) 'Fourier constants'
480             do k=1,nterm(i,j)
481               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
482             enddo
483             write (iout,*) 'Lorenz constants'
484             do k=1,nlor(i,j)
485               write (iout,'(3(1pe15.5))') 
486      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
487             enddo
488           enddo
489         enddo
490       endif
491 C
492 C 6/23/01 Read parameters for double torsionals
493 C
494       do i=1,ntortyp
495         do j=1,ntortyp
496           do k=1,ntortyp
497             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
498             if (t1.ne.onelett(i) .or. t2.ne.onelett(j) 
499      &        .or. t3.ne.onelett(k)) then
500               write (iout,*) "Error in double torsional parameter file",
501      &         i,j,k,t1,t2,t3
502 #ifdef MPI
503               call MPI_Finalize(Ierror)
504 #endif
505                stop "Error in double torsional parameter file"
506             endif
507             read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
508      &         ntermd_2(i,j,k)
509             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
510      &         ntermd_1(i,j,k))
511             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
512      &         ntermd_1(i,j,k))
513             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
514      &         ntermd_1(i,j,k))
515             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
516      &         ntermd_1(i,j,k))
517             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
518      &         v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
519      &         m=1,l-1),l=1,ntermd_2(i,j,k))
520           enddo
521         enddo
522       enddo
523       if (lprint) then
524       write (iout,*) 
525       write (iout,*) 'Constants for double torsionals'
526       do i=1,ntortyp
527         do j=1,ntortyp 
528           do k=1,ntortyp
529             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
530      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
531             write (iout,*)
532             write (iout,*) 'Single angles:'
533             do l=1,ntermd_1(i,j,k)
534               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
535      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
536      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
537             enddo
538             write (iout,*)
539             write (iout,*) 'Pairs of angles:'
540             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
541             do l=1,ntermd_2(i,j,k)
542               write (iout,'(i5,20f10.5)') 
543      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
544             enddo
545             write (iout,*)
546             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
547             do l=1,ntermd_2(i,j,k)
548               write (iout,'(i5,20f10.5)') 
549      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
550             enddo
551             write (iout,*)
552           enddo
553         enddo
554       enddo
555       endif
556 #endif
557 C Read of Side-chain backbone correlation parameters
558 C Modified 11 May 2012 by Adasko
559 CCC
560 C
561 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
562 C
563       read (isccor,*,end=119,err=119) nsccortyp
564 #ifdef SCCORPDB
565       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
566       do i=-ntyp,-1
567         isccortyp(i)=-isccortyp(-i)
568       enddo
569       iscprol=isccortyp(20)
570 c      write (iout,*) 'ntortyp',ntortyp
571       maxinter=3
572 cc maxinter is maximum interaction sites
573       do l=1,maxinter
574       do i=1,nsccortyp
575         do j=1,nsccortyp
576           read (isccor,*,end=119,err=119)
577      &nterm_sccor(i,j),nlor_sccor(i,j)
578           v0ijsccor=0.0d0
579           v0ijsccor1=0.0d0
580           v0ijsccor2=0.0d0
581           v0ijsccor3=0.0d0
582           si=-1.0d0
583           nterm_sccor(-i,j)=nterm_sccor(i,j)
584           nterm_sccor(-i,-j)=nterm_sccor(i,j)
585           nterm_sccor(i,-j)=nterm_sccor(i,j)
586           do k=1,nterm_sccor(i,j)
587             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
588      &    ,v2sccor(k,l,i,j)
589             if (j.eq.iscprol) then
590              if (i.eq.isccortyp(10)) then
591              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
592              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
593              else
594              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
595      &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
596              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
597      &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
598              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
599              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
600              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
601              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
602              endif
603             else
604              if (i.eq.isccortyp(10)) then
605              v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
606              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
607              else
608                if (j.eq.isccortyp(10)) then
609              v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
610              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
611                else
612              v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
613              v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
614              v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
615              v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
616              v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
617              v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
618                 endif
619                endif
620             endif
621             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
622             v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
623             v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
624             v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
625             si=-si
626           enddo
627           do k=1,nlor_sccor(i,j)
628             read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
629      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
630             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
631      &(1+vlor3sccor(k,i,j)**2)
632           enddo
633           v0sccor(l,i,j)=v0ijsccor
634           v0sccor(l,-i,j)=v0ijsccor1
635           v0sccor(l,i,-j)=v0ijsccor2
636           v0sccor(l,-i,-j)=v0ijsccor3         
637         enddo
638       enddo
639       enddo
640       close (isccor)
641 #else
642       read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
643 c      write (iout,*) 'ntortyp',ntortyp
644       maxinter=3
645 cc maxinter is maximum interaction sites
646       do l=1,maxinter    
647       do i=1,nsccortyp
648         do j=1,nsccortyp
649           read (isccor,*,end=113,err=113)
650      & nterm_sccor(i,j),nlor_sccor(i,j)
651           v0ijsccor=0.0d0
652           si=-1.0d0
653   
654           do k=1,nterm_sccor(i,j)
655             read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
656      &    ,v2sccor(k,l,i,j) 
657             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
658             si=-si
659           enddo
660           do k=1,nlor_sccor(i,j)
661             read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
662      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j) 
663             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
664      &(1+vlor3sccor(k,i,j)**2)
665           enddo
666           v0sccor(l,i,j)=v0ijsccor
667         enddo
668       enddo
669       enddo
670       close (isccor)
671       
672 #endif      
673       if (lprint) then
674         write (iout,'(/a/)') 'Torsional constants:'
675         do i=1,nsccortyp
676           do j=1,nsccortyp
677             write (iout,*) 'ityp',i,' jtyp',j
678             write (iout,*) 'Fourier constants'
679             do k=1,nterm_sccor(i,j)
680       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
681             enddo
682             write (iout,*) 'Lorenz constants'
683             do k=1,nlor_sccor(i,j)
684               write (iout,'(3(1pe15.5))') 
685      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
686             enddo
687           enddo
688         enddo
689       endif
690
691 C
692 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
693 C         interaction energy of the Gly, Ala, and Pro prototypes.
694 C
695       if (lprint) then
696         write (iout,*)
697         write (iout,*) "Coefficients of the cumulants"
698       endif
699       read (ifourier,*) nloctyp
700       do i=1,nloctyp
701         read (ifourier,*,end=115,err=115)
702         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
703         if (lprint) then
704         write (iout,*) 'Type',i
705         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
706         endif
707         B1(1,i)  = b(3)
708         B1(2,i)  = b(5)
709 c        b1(1,i)=0.0d0
710 c        b1(2,i)=0.0d0
711         B1tilde(1,i) = b(3)
712         B1tilde(2,i) =-b(5) 
713 c        b1tilde(1,i)=0.0d0
714 c        b1tilde(2,i)=0.0d0
715         B2(1,i)  = b(2)
716         B2(2,i)  = b(4)
717 c        b2(1,i)=0.0d0
718 c        b2(2,i)=0.0d0
719         CC(1,1,i)= b(7)
720         CC(2,2,i)=-b(7)
721         CC(2,1,i)= b(9)
722         CC(1,2,i)= b(9)
723 c        CC(1,1,i)=0.0d0
724 c        CC(2,2,i)=0.0d0
725 c        CC(2,1,i)=0.0d0
726 c        CC(1,2,i)=0.0d0
727         Ctilde(1,1,i)=b(7)
728         Ctilde(1,2,i)=b(9)
729         Ctilde(2,1,i)=-b(9)
730         Ctilde(2,2,i)=b(7)
731 c        Ctilde(1,1,i)=0.0d0
732 c        Ctilde(1,2,i)=0.0d0
733 c        Ctilde(2,1,i)=0.0d0
734 c        Ctilde(2,2,i)=0.0d0
735         DD(1,1,i)= b(6)
736         DD(2,2,i)=-b(6)
737         DD(2,1,i)= b(8)
738         DD(1,2,i)= b(8)
739 c        DD(1,1,i)=0.0d0
740 c        DD(2,2,i)=0.0d0
741 c        DD(2,1,i)=0.0d0
742 c        DD(1,2,i)=0.0d0
743         Dtilde(1,1,i)=b(6)
744         Dtilde(1,2,i)=b(8)
745         Dtilde(2,1,i)=-b(8)
746         Dtilde(2,2,i)=b(6)
747 c        Dtilde(1,1,i)=0.0d0
748 c        Dtilde(1,2,i)=0.0d0
749 c        Dtilde(2,1,i)=0.0d0
750 c        Dtilde(2,2,i)=0.0d0
751         EE(1,1,i)= b(10)+b(11)
752         EE(2,2,i)=-b(10)+b(11)
753         EE(2,1,i)= b(12)-b(13)
754         EE(1,2,i)= b(12)+b(13)
755 c        ee(1,1,i)=1.0d0
756 c        ee(2,2,i)=1.0d0
757 c        ee(2,1,i)=0.0d0
758 c        ee(1,2,i)=0.0d0
759 c        ee(2,1,i)=ee(1,2,i)
760       enddo
761 C      write(iout,*) "parm", B1(1,nloctyp+1),B1(2,nloctyp+1)
762       if (lprint) then
763       do i=1,nloctyp+1
764         write (iout,*) 'Type',i
765         write (iout,*) 'B1'
766         write(iout,*) B1(1,i),B1(2,i)
767         write (iout,*) 'B2'
768         write(iout,*) B2(1,i),B2(2,i)
769         write (iout,*) 'CC'
770         do j=1,2
771           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
772         enddo
773         write(iout,*) 'DD'
774         do j=1,2
775           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
776         enddo
777         write(iout,*) 'EE'
778         do j=1,2
779           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
780         enddo
781       enddo
782       endif
783
784 C Read electrostatic-interaction parameters
785 C
786       if (lprint) then
787         write (iout,*)
788         write (iout,'(/a)') 'Electrostatic interaction constants:'
789         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') 
790      &            'IT','JT','APP','BPP','AEL6','AEL3'
791       endif
792       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
793       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
794       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
795       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
796       close (ielep)
797       do i=1,2
798         do j=1,2
799         rri=rpp(i,j)**6
800         app (i,j)=epp(i,j)*rri*rri 
801         bpp (i,j)=-2.0D0*epp(i,j)*rri
802         ael6(i,j)=elpp6(i,j)*4.2D0**6
803         ael3(i,j)=elpp3(i,j)*4.2D0**3
804         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
805      &                    ael6(i,j),ael3(i,j)
806         enddo
807       enddo
808 C
809 C Read side-chain interaction parameters.
810 C
811       read (isidep,*,end=117,err=117) ipot,expon
812       if (ipot.lt.1 .or. ipot.gt.5) then
813         write (iout,'(2a)') 'Error while reading SC interaction',
814      &               'potential file - unknown potential type.'
815 #ifdef MPI
816         call MPI_Finalize(Ierror)
817 #endif
818         stop
819       endif
820       expon2=expon/2
821       if(me.eq.king)
822      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
823      & ', exponents are ',expon,2*expon 
824       goto (10,20,30,30,40) ipot
825 C----------------------- LJ potential ---------------------------------
826    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
827      &   (sigma0(i),i=1,ntyp)
828       if (lprint) then
829         write (iout,'(/a/)') 'Parameters of the LJ potential:'
830         write (iout,'(a/)') 'The epsilon array:'
831         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
832         write (iout,'(/a)') 'One-body parameters:'
833         write (iout,'(a,4x,a)') 'residue','sigma'
834         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
835       endif
836       goto 50
837 C----------------------- LJK potential --------------------------------
838    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
839      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
840       if (lprint) then
841         write (iout,'(/a/)') 'Parameters of the LJK potential:'
842         write (iout,'(a/)') 'The epsilon array:'
843         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
844         write (iout,'(/a)') 'One-body parameters:'
845         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
846         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
847      &        i=1,ntyp)
848       endif
849       goto 50
850 C---------------------- GB or BP potential -----------------------------
851    30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
852      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
853      &  (alp(i),i=1,ntyp)
854 C For the GB potential convert sigma'**2 into chi'
855       if (ipot.eq.4) then
856         do i=1,ntyp
857           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
858         enddo
859       endif
860       if (lprint) then
861         write (iout,'(/a/)') 'Parameters of the BP potential:'
862         write (iout,'(a/)') 'The epsilon array:'
863         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
864         write (iout,'(/a)') 'One-body parameters:'
865         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
866      &       '    chip  ','    alph  '
867         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
868      &                     chip(i),alp(i),i=1,ntyp)
869       endif
870       goto 50
871 C--------------------- GBV potential -----------------------------------
872    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
873      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
874      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
875       if (lprint) then
876         write (iout,'(/a/)') 'Parameters of the GBV potential:'
877         write (iout,'(a/)') 'The epsilon array:'
878         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
879         write (iout,'(/a)') 'One-body parameters:'
880         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
881      &      's||/s_|_^2','    chip  ','    alph  '
882         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
883      &           sigii(i),chip(i),alp(i),i=1,ntyp)
884       endif
885    50 continue
886       close (isidep)
887 C-----------------------------------------------------------------------
888 C Calculate the "working" parameters of SC interactions.
889       do i=2,ntyp
890         do j=1,i-1
891           eps(i,j)=eps(j,i)
892         enddo
893       enddo
894       do i=1,ntyp
895         do j=i,ntyp
896           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
897           sigma(j,i)=sigma(i,j)
898           rs0(i,j)=dwa16*sigma(i,j)
899           rs0(j,i)=rs0(i,j)
900         enddo
901       enddo
902       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') 
903      & 'Working parameters of the SC interactions:',
904      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
905      & '  chi1   ','   chi2   ' 
906       do i=1,ntyp
907         do j=i,ntyp
908           epsij=eps(i,j)
909           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
910             rrij=sigma(i,j)
911           else
912             rrij=rr0(i)+rr0(j)
913           endif
914           r0(i,j)=rrij
915           r0(j,i)=rrij
916           rrij=rrij**expon
917           epsij=eps(i,j)
918           sigeps=dsign(1.0D0,epsij)
919           epsij=dabs(epsij)
920           aa(i,j)=epsij*rrij*rrij
921           bb(i,j)=-sigeps*epsij*rrij
922           aa(j,i)=aa(i,j)
923           bb(j,i)=bb(i,j)
924           if (ipot.gt.2) then
925             sigt1sq=sigma0(i)**2
926             sigt2sq=sigma0(j)**2
927             sigii1=sigii(i)
928             sigii2=sigii(j)
929             ratsig1=sigt2sq/sigt1sq
930             ratsig2=1.0D0/ratsig1
931             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
932             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
933             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
934           else
935             rsum_max=sigma(i,j)
936           endif
937 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
938             sigmaii(i,j)=rsum_max
939             sigmaii(j,i)=rsum_max 
940 c         else
941 c           sigmaii(i,j)=r0(i,j)
942 c           sigmaii(j,i)=r0(i,j)
943 c         endif
944 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
945           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
946             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
947             augm(i,j)=epsij*r_augm**(2*expon)
948 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
949             augm(j,i)=augm(i,j)
950           else
951             augm(i,j)=0.0D0
952             augm(j,i)=0.0D0
953           endif
954           if (lprint) then
955             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
956      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
957      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
958           endif
959         enddo
960       enddo
961 #ifdef OLDSCP
962 C
963 C Define the SC-p interaction constants (hard-coded; old style)
964 C
965       do i=1,20
966 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
967 C helix formation)
968 c       aad(i,1)=0.3D0*4.0D0**12
969 C Following line for constants currently implemented
970 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
971         aad(i,1)=1.5D0*4.0D0**12
972 c       aad(i,1)=0.17D0*5.6D0**12
973         aad(i,2)=aad(i,1)
974 C "Soft" SC-p repulsion
975         bad(i,1)=0.0D0
976 C Following line for constants currently implemented
977 c       aad(i,1)=0.3D0*4.0D0**6
978 C "Hard" SC-p repulsion
979         bad(i,1)=3.0D0*4.0D0**6
980 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
981         bad(i,2)=bad(i,1)
982 c       aad(i,1)=0.0D0
983 c       aad(i,2)=0.0D0
984 c       bad(i,1)=1228.8D0
985 c       bad(i,2)=1228.8D0
986       enddo
987 #else
988 C
989 C 8/9/01 Read the SC-p interaction constants from file
990 C
991       do i=1,ntyp
992         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
993       enddo
994       do i=1,ntyp
995         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
996         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
997         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
998         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
999       enddo
1000
1001       if (lprint) then
1002         write (iout,*) "Parameters of SC-p interactions:"
1003         do i=1,20
1004           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
1005      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
1006         enddo
1007       endif
1008 #endif
1009 C
1010 C Define the constants of the disulfide bridge
1011 C
1012       ebr=-5.50D0
1013 c
1014 c Old arbitrary potential - commented out.
1015 c
1016 c      dbr= 4.20D0
1017 c      fbr= 3.30D0
1018 c
1019 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1020 c energy surface of diethyl disulfide.
1021 c A. Liwo and U. Kozlowska, 11/24/03
1022 c
1023       D0CM = 3.78d0
1024       AKCM = 15.1d0
1025       AKTH = 11.0d0
1026       AKCT = 12.0d0
1027       V1SS =-1.08d0
1028       V2SS = 7.61d0
1029       V3SS = 13.7d0
1030 c      akcm=0.0d0
1031 c      akth=0.0d0
1032 c      akct=0.0d0
1033 c      v1ss=0.0d0
1034 c      v2ss=0.0d0
1035 c      v3ss=0.0d0
1036       
1037       if(me.eq.king) then
1038       write (iout,'(/a)') "Disulfide bridge parameters:"
1039       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1040       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1041       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1042       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1043      &  ' v3ss:',v3ss
1044       endif
1045       return
1046   111 write (iout,*) "Error reading bending energy parameters."
1047       goto 999
1048   112 write (iout,*) "Error reading rotamer energy parameters."
1049       goto 999
1050   113 write (iout,*) "Error reading torsional energy parameters."
1051       goto 999
1052   114 write (iout,*) "Error reading double torsional energy parameters."
1053       goto 999
1054   115 write (iout,*) 
1055      &  "Error reading cumulant (multibody energy) parameters."
1056       goto 999
1057   116 write (iout,*) "Error reading electrostatic energy parameters."
1058       goto 999
1059   117 write (iout,*) "Error reading side chain interaction parameters."
1060       goto 999
1061   118 write (iout,*) "Error reading SCp interaction parameters."
1062       goto 999
1063   119 write (iout,*) "Error reading SCCOR parameters"
1064   999 continue
1065 #ifdef MPI
1066       call MPI_Finalize(Ierror)
1067 #endif
1068       stop
1069       return
1070       end
1071
1072
1073       subroutine getenv_loc(var, val)
1074       character(*) var, val
1075
1076 #ifdef WINIFL
1077       character(2000) line
1078       external ilen
1079
1080       open (196,file='env',status='old',readonly,shared)
1081       iread=0
1082 c      write(*,*)'looking for ',var
1083 10    read(196,*,err=11,end=11)line
1084       iread=index(line,var)
1085 c      write(*,*)iread,' ',var,' ',line
1086       if (iread.eq.0) go to 10 
1087 c      write(*,*)'---> ',line
1088 11    continue
1089       if(iread.eq.0) then
1090 c       write(*,*)'CHUJ'
1091        val=''
1092       else
1093        iread=iread+ilen(var)+1
1094        read (line(iread:),*,err=12,end=12) val
1095 c       write(*,*)'OK: ',var,' = ',val
1096       endif
1097       close(196)
1098       return
1099 12    val=''
1100       close(196)
1101 #elif (defined CRAY)
1102       integer lennam,lenval,ierror
1103 c
1104 c        getenv using a POSIX call, useful on the T3D
1105 c        Sept 1996, comment out error check on advice of H. Pritchard
1106 c
1107       lennam = len(var)
1108       if(lennam.le.0) stop '--error calling getenv--'
1109       call pxfgetenv(var,lennam,val,lenval,ierror)
1110 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1111 #else
1112       call getenv(var,val)
1113 #endif
1114
1115       return
1116       end