ctest the first static disulfide test ene/min cart
[unres.git] / source / unres / src_MIN / 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 c      write (iout,*) 'ntortyp',ntortyp
454       do i=1,ntortyp
455         do j=1,ntortyp
456           read (itorp,*,end=113,err=113) nterm(i,j),nlor(i,j)
457           v0ij=0.0d0
458           si=-1.0d0
459           do k=1,nterm(i,j)
460             read (itorp,*,end=113,err=113) kk,v1(k,i,j),v2(k,i,j) 
461             v0ij=v0ij+si*v1(k,i,j)
462             si=-si
463           enddo
464           do k=1,nlor(i,j)
465             read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),
466      &        vlor2(k,i,j),vlor3(k,i,j) 
467             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
468           enddo
469           v0(i,j)=v0ij
470         enddo
471       enddo
472       close (itorp)
473       if (lprint) then
474         write (iout,'(/a/)') 'Torsional constants:'
475         do i=1,ntortyp
476           do j=1,ntortyp
477             write (iout,*) 'ityp',i,' jtyp',j
478             write (iout,*) 'Fourier constants'
479             do k=1,nterm(i,j)
480               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
481             enddo
482             write (iout,*) 'Lorenz constants'
483             do k=1,nlor(i,j)
484               write (iout,'(3(1pe15.5))') 
485      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
486             enddo
487           enddo
488         enddo
489       endif
490 C
491 C 6/23/01 Read parameters for double torsionals
492 C
493       do i=1,ntortyp
494         do j=1,ntortyp
495           do k=1,ntortyp
496             read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
497             if (t1.ne.onelett(i) .or. t2.ne.onelett(j) 
498      &        .or. t3.ne.onelett(k)) then
499               write (iout,*) "Error in double torsional parameter file",
500      &         i,j,k,t1,t2,t3
501 #ifdef MPI
502               call MPI_Finalize(Ierror)
503 #endif
504                stop "Error in double torsional parameter file"
505             endif
506             read (itordp,*,end=114,err=114) ntermd_1(i,j,k),
507      &         ntermd_2(i,j,k)
508             read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k),l=1,
509      &         ntermd_1(i,j,k))
510             read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k),l=1,
511      &         ntermd_1(i,j,k))
512             read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k),l=1,
513      &         ntermd_1(i,j,k))
514             read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k),l=1,
515      &         ntermd_1(i,j,k))
516             read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k),
517      &         v2c(m,l,i,j,k),v2s(l,m,i,j,k),v2s(m,l,i,j,k),
518      &         m=1,l-1),l=1,ntermd_2(i,j,k))
519           enddo
520         enddo
521       enddo
522       if (lprint) then
523       write (iout,*) 
524       write (iout,*) 'Constants for double torsionals'
525       do i=1,ntortyp
526         do j=1,ntortyp 
527           do k=1,ntortyp
528             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
529      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
530             write (iout,*)
531             write (iout,*) 'Single angles:'
532             do l=1,ntermd_1(i,j,k)
533               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
534      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
535      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
536             enddo
537             write (iout,*)
538             write (iout,*) 'Pairs of angles:'
539             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
540             do l=1,ntermd_2(i,j,k)
541               write (iout,'(i5,20f10.5)') 
542      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
543             enddo
544             write (iout,*)
545             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
546             do l=1,ntermd_2(i,j,k)
547               write (iout,'(i5,20f10.5)') 
548      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
549             enddo
550             write (iout,*)
551           enddo
552         enddo
553       enddo
554       endif
555 #endif
556 C
557 C 5/21/07 (AL) Read coefficients of the backbone-local sidechain-local
558 C         correlation energies.
559 C
560       read (isccor,*,end=119,err=119) nsccortyp
561 #ifdef SCCORPDB
562       read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
563 C For D-aminoacid uncomment
564 C      do i=-ntyp,-1
565 C        isccortyp(i)=-isccortyp(-i)
566 C      enddo
567       iscprol=isccortyp(20)
568 c      write (iout,*) 'ntortyp',ntortyp
569       maxinter=3
570 cc maxinter is maximum interaction sites
571       do l=1,maxinter
572       do i=1,nsccortyp
573         do j=1,nsccortyp
574           read (isccor,*,end=119,err=119)
575      &nterm_sccor(i,j),nlor_sccor(i,j)
576           v0ijsccor=0.0d0
577           si=-1.0d0
578 C          nterm_sccor(-i,j)=nterm_sccor(i,j)
579 C          nterm_sccor(-i,-j)=nterm_sccor(i,j)
580 C          nterm_sccor(i,-j)=nterm_sccor(i,j)
581           do k=1,nterm_sccor(i,j)
582             read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j)
583      &    ,v2sccor(k,l,i,j)
584 C            if (j.eq.iscprol) then
585 C             if (i.eq.isccortyp(10)) then
586 C             v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
587 C             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
588 C             else
589 C             v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0
590 C     &                        +v2sccor(k,l,i,j)*dsqrt(0.75d0)
591 C             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0
592 C     &                        +v1sccor(k,l,i,j)*dsqrt(0.75d0)
593 C             v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
594 C             v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
595 C             v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
596 C             v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
597 C             endif
598 CC           else
599 C             if (i.eq.isccortyp(10)) then
600 C             v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
601 C             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
602 C             else
603 C               if (j.eq.isccortyp(10)) then
604 C             v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
605 C             v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
606 C               else
607 C             v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
608 C             v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
609 C             v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
610 C             v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
611 C             v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
612 C             v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
613 C                endif
614 C               endif
615 C            endif
616             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
617             si=-si
618           enddo
619           do k=1,nlor_sccor(i,j)
620             read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
621      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
622             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
623      &(1+vlor3sccor(k,i,j)**2)
624           enddo
625           v0sccor(i,j)=v0ijsccor
626         enddo
627       enddo
628       enddo
629       close (isccor)
630 #else
631       read (isccor,*,end=113,err=113) (isccortyp(i),i=1,ntyp)
632 c      write (iout,*) 'ntortyp',ntortyp
633       maxinter=3
634 cc maxinter is maximum interaction sites
635       do l=1,maxinter
636       do i=1,nsccortyp
637         do j=1,nsccortyp
638           read (isccor,*,end=113,err=113)
639      & nterm_sccor(i,j),nlor_sccor(i,j)
640           v0ijsccor=0.0d0
641           si=-1.0d0
642
643           do k=1,nterm_sccor(i,j)
644             read (isccor,*,end=113,err=113) kk,v1sccor(k,l,i,j)
645      &    ,v2sccor(k,l,i,j)
646             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
647             si=-si
648           enddo
649           do k=1,nlor_sccor(i,j)
650             read (isccor,*,end=113,err=113) kk,vlor1sccor(k,i,j),
651      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
652             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
653      &(1+vlor3sccor(k,i,j)**2)
654           enddo
655           v0sccor(i,j)=v0ijsccor
656         enddo
657       enddo
658       enddo
659       close (isccor)
660
661 #endif      
662       if (lprint) then
663        write (iout,'(/a/)') 'Torsional constants of SCCORR:'
664 C       do i=1,20
665 C         do j=1,20
666 C        write (iout,'(/a/)') 'Torsional constants:'
667         do i=1,nsccortyp
668           do j=1,nsccortyp
669             write (iout,*) 'ityp',i,' jtyp',j
670 C            do k=1,nterm_sccor
671 C             write (iout,'(2(1pe15.5))') v1sccor(k,i,j),v2sccor(k,i,j)
672 C            write (iout,*) 'Fourier constants'
673             do k=1,nterm_sccor(i,j)
674       write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
675             enddo
676             write (iout,*) 'Lorenz constants'
677             do k=1,nlor_sccor(i,j)
678               write (iout,'(3(1pe15.5))')
679      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
680             enddo
681           enddo
682         enddo
683       endif
684 C
685 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
686 C         interaction energy of the Gly, Ala, and Pro prototypes.
687 C
688       if (lprint) then
689         write (iout,*)
690         write (iout,*) "Coefficients of the cumulants"
691       endif
692       read (ifourier,*) nloctyp
693       do i=1,nloctyp
694         read (ifourier,*,end=115,err=115)
695         read (ifourier,*,end=115,err=115) (b(ii),ii=1,13)
696         if (lprint) then
697         write (iout,*) 'Type',i
698         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii),ii=1,13)
699         endif
700         B1(1,i)  = b(3)
701         B1(2,i)  = b(5)
702 c        b1(1,i)=0.0d0
703 c        b1(2,i)=0.0d0
704         B1tilde(1,i) = b(3)
705         B1tilde(2,i) =-b(5) 
706 c        b1tilde(1,i)=0.0d0
707 c        b1tilde(2,i)=0.0d0
708         B2(1,i)  = b(2)
709         B2(2,i)  = b(4)
710 c        b2(1,i)=0.0d0
711 c        b2(2,i)=0.0d0
712         CC(1,1,i)= b(7)
713         CC(2,2,i)=-b(7)
714         CC(2,1,i)= b(9)
715         CC(1,2,i)= b(9)
716 c        CC(1,1,i)=0.0d0
717 c        CC(2,2,i)=0.0d0
718 c        CC(2,1,i)=0.0d0
719 c        CC(1,2,i)=0.0d0
720         Ctilde(1,1,i)=b(7)
721         Ctilde(1,2,i)=b(9)
722         Ctilde(2,1,i)=-b(9)
723         Ctilde(2,2,i)=b(7)
724 c        Ctilde(1,1,i)=0.0d0
725 c        Ctilde(1,2,i)=0.0d0
726 c        Ctilde(2,1,i)=0.0d0
727 c        Ctilde(2,2,i)=0.0d0
728         DD(1,1,i)= b(6)
729         DD(2,2,i)=-b(6)
730         DD(2,1,i)= b(8)
731         DD(1,2,i)= b(8)
732 c        DD(1,1,i)=0.0d0
733 c        DD(2,2,i)=0.0d0
734 c        DD(2,1,i)=0.0d0
735 c        DD(1,2,i)=0.0d0
736         Dtilde(1,1,i)=b(6)
737         Dtilde(1,2,i)=b(8)
738         Dtilde(2,1,i)=-b(8)
739         Dtilde(2,2,i)=b(6)
740 c        Dtilde(1,1,i)=0.0d0
741 c        Dtilde(1,2,i)=0.0d0
742 c        Dtilde(2,1,i)=0.0d0
743 c        Dtilde(2,2,i)=0.0d0
744         EE(1,1,i)= b(10)+b(11)
745         EE(2,2,i)=-b(10)+b(11)
746         EE(2,1,i)= b(12)-b(13)
747         EE(1,2,i)= b(12)+b(13)
748 c        ee(1,1,i)=1.0d0
749 c        ee(2,2,i)=1.0d0
750 c        ee(2,1,i)=0.0d0
751 c        ee(1,2,i)=0.0d0
752 c        ee(2,1,i)=ee(1,2,i)
753       enddo
754       if (lprint) then
755       do i=1,nloctyp
756         write (iout,*) 'Type',i
757         write (iout,*) 'B1'
758         write(iout,*) B1(1,i),B1(2,i)
759         write (iout,*) 'B2'
760         write(iout,*) B2(1,i),B2(2,i)
761         write (iout,*) 'CC'
762         do j=1,2
763           write (iout,'(2f10.5)') CC(j,1,i),CC(j,2,i)
764         enddo
765         write(iout,*) 'DD'
766         do j=1,2
767           write (iout,'(2f10.5)') DD(j,1,i),DD(j,2,i)
768         enddo
769         write(iout,*) 'EE'
770         do j=1,2
771           write (iout,'(2f10.5)') EE(j,1,i),EE(j,2,i)
772         enddo
773       enddo
774       endif
775
776 C Read electrostatic-interaction parameters
777 C
778       if (lprint) then
779         write (iout,*)
780         write (iout,'(/a)') 'Electrostatic interaction constants:'
781         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') 
782      &            'IT','JT','APP','BPP','AEL6','AEL3'
783       endif
784       read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
785       read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
786       read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
787       read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
788       close (ielep)
789       do i=1,2
790         do j=1,2
791         rri=rpp(i,j)**6
792         app (i,j)=epp(i,j)*rri*rri 
793         bpp (i,j)=-2.0D0*epp(i,j)*rri
794         ael6(i,j)=elpp6(i,j)*4.2D0**6
795         ael3(i,j)=elpp3(i,j)*4.2D0**3
796         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
797      &                    ael6(i,j),ael3(i,j)
798         enddo
799       enddo
800 C
801 C Read side-chain interaction parameters.
802 C
803       read (isidep,*,end=117,err=117) ipot,expon
804       if (ipot.lt.1 .or. ipot.gt.5) then
805         write (iout,'(2a)') 'Error while reading SC interaction',
806      &               'potential file - unknown potential type.'
807 #ifdef MPI
808         call MPI_Finalize(Ierror)
809 #endif
810         stop
811       endif
812       expon2=expon/2
813       if(me.eq.king)
814      & write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
815      & ', exponents are ',expon,2*expon 
816       goto (10,20,30,30,40) ipot
817 C----------------------- LJ potential ---------------------------------
818    10 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
819      &   (sigma0(i),i=1,ntyp)
820       if (lprint) then
821         write (iout,'(/a/)') 'Parameters of the LJ potential:'
822         write (iout,'(a/)') 'The epsilon array:'
823         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
824         write (iout,'(/a)') 'One-body parameters:'
825         write (iout,'(a,4x,a)') 'residue','sigma'
826         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
827       endif
828       goto 50
829 C----------------------- LJK potential --------------------------------
830    20 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
831      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
832       if (lprint) then
833         write (iout,'(/a/)') 'Parameters of the LJK potential:'
834         write (iout,'(a/)') 'The epsilon array:'
835         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
836         write (iout,'(/a)') 'One-body parameters:'
837         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
838         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
839      &        i=1,ntyp)
840       endif
841       goto 50
842 C---------------------- GB or BP potential -----------------------------
843    30 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
844      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip(i),i=1,ntyp),
845      &  (alp(i),i=1,ntyp)
846 C For the GB potential convert sigma'**2 into chi'
847       if (ipot.eq.4) then
848         do i=1,ntyp
849           chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
850         enddo
851       endif
852       if (lprint) then
853         write (iout,'(/a/)') 'Parameters of the BP potential:'
854         write (iout,'(a/)') 'The epsilon array:'
855         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
856         write (iout,'(/a)') 'One-body parameters:'
857         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
858      &       '    chip  ','    alph  '
859         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
860      &                     chip(i),alp(i),i=1,ntyp)
861       endif
862       goto 50
863 C--------------------- GBV potential -----------------------------------
864    40 read (isidep,*,end=116,err=116)((eps(i,j),j=i,ntyp),i=1,ntyp),
865      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
866      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
867       if (lprint) then
868         write (iout,'(/a/)') 'Parameters of the GBV potential:'
869         write (iout,'(a/)') 'The epsilon array:'
870         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
871         write (iout,'(/a)') 'One-body parameters:'
872         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
873      &      's||/s_|_^2','    chip  ','    alph  '
874         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
875      &           sigii(i),chip(i),alp(i),i=1,ntyp)
876       endif
877    50 continue
878       close (isidep)
879 C-----------------------------------------------------------------------
880 C Calculate the "working" parameters of SC interactions.
881       do i=2,ntyp
882         do j=1,i-1
883           eps(i,j)=eps(j,i)
884         enddo
885       enddo
886       do i=1,ntyp
887         do j=i,ntyp
888           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
889           sigma(j,i)=sigma(i,j)
890           rs0(i,j)=dwa16*sigma(i,j)
891           rs0(j,i)=rs0(i,j)
892         enddo
893       enddo
894       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') 
895      & 'Working parameters of the SC interactions:',
896      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
897      & '  chi1   ','   chi2   ' 
898       do i=1,ntyp
899         do j=i,ntyp
900           epsij=eps(i,j)
901           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
902             rrij=sigma(i,j)
903           else
904             rrij=rr0(i)+rr0(j)
905           endif
906           r0(i,j)=rrij
907           r0(j,i)=rrij
908           rrij=rrij**expon
909           epsij=eps(i,j)
910           sigeps=dsign(1.0D0,epsij)
911           epsij=dabs(epsij)
912           aa(i,j)=epsij*rrij*rrij
913           bb(i,j)=-sigeps*epsij*rrij
914           aa(j,i)=aa(i,j)
915           bb(j,i)=bb(i,j)
916           if (ipot.gt.2) then
917             sigt1sq=sigma0(i)**2
918             sigt2sq=sigma0(j)**2
919             sigii1=sigii(i)
920             sigii2=sigii(j)
921             ratsig1=sigt2sq/sigt1sq
922             ratsig2=1.0D0/ratsig1
923             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
924             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
925             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
926           else
927             rsum_max=sigma(i,j)
928           endif
929 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
930             sigmaii(i,j)=rsum_max
931             sigmaii(j,i)=rsum_max 
932 c         else
933 c           sigmaii(i,j)=r0(i,j)
934 c           sigmaii(j,i)=r0(i,j)
935 c         endif
936 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
937           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
938             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
939             augm(i,j)=epsij*r_augm**(2*expon)
940 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
941             augm(j,i)=augm(i,j)
942           else
943             augm(i,j)=0.0D0
944             augm(j,i)=0.0D0
945           endif
946           if (lprint) then
947             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
948      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
949      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
950           endif
951         enddo
952       enddo
953 #ifdef OLDSCP
954 C
955 C Define the SC-p interaction constants (hard-coded; old style)
956 C
957       do i=1,20
958 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
959 C helix formation)
960 c       aad(i,1)=0.3D0*4.0D0**12
961 C Following line for constants currently implemented
962 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
963         aad(i,1)=1.5D0*4.0D0**12
964 c       aad(i,1)=0.17D0*5.6D0**12
965         aad(i,2)=aad(i,1)
966 C "Soft" SC-p repulsion
967         bad(i,1)=0.0D0
968 C Following line for constants currently implemented
969 c       aad(i,1)=0.3D0*4.0D0**6
970 C "Hard" SC-p repulsion
971         bad(i,1)=3.0D0*4.0D0**6
972 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
973         bad(i,2)=bad(i,1)
974 c       aad(i,1)=0.0D0
975 c       aad(i,2)=0.0D0
976 c       bad(i,1)=1228.8D0
977 c       bad(i,2)=1228.8D0
978       enddo
979 #else
980 C
981 C 8/9/01 Read the SC-p interaction constants from file
982 C
983       do i=1,ntyp
984         read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
985       enddo
986       do i=1,ntyp
987         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
988         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
989         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
990         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
991       enddo
992
993       if (lprint) then
994         write (iout,*) "Parameters of SC-p interactions:"
995         do i=1,20
996           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
997      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
998         enddo
999       endif
1000 #endif
1001 C
1002 C Define the constants of the disulfide bridge
1003 C
1004       ebr=-5.50D0
1005 c
1006 c Old arbitrary potential - commented out.
1007 c
1008 c      dbr= 4.20D0
1009 c      fbr= 3.30D0
1010 c
1011 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
1012 c energy surface of diethyl disulfide.
1013 c A. Liwo and U. Kozlowska, 11/24/03
1014 c
1015       D0CM = 3.78d0
1016       AKCM = 15.1d0
1017       AKTH = 11.0d0
1018       AKCT = 12.0d0
1019       V1SS =-1.08d0
1020       V2SS = 7.61d0
1021       V3SS = 13.7d0
1022 c      akcm=0.0d0
1023 c      akth=0.0d0
1024 c      akct=0.0d0
1025 c      v1ss=0.0d0
1026 c      v2ss=0.0d0
1027 c      v3ss=0.0d0
1028       
1029       if(me.eq.king) then
1030       write (iout,'(/a)') "Disulfide bridge parameters:"
1031       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
1032       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
1033       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
1034       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
1035      &  ' v3ss:',v3ss
1036       endif
1037       return
1038   111 write (iout,*) "Error reading bending energy parameters."
1039       goto 999
1040   112 write (iout,*) "Error reading rotamer energy parameters."
1041       goto 999
1042   113 write (iout,*) "Error reading torsional energy parameters."
1043       goto 999
1044   114 write (iout,*) "Error reading double torsional energy parameters."
1045       goto 999
1046   115 write (iout,*) 
1047      &  "Error reading cumulant (multibody energy) parameters."
1048       goto 999
1049   116 write (iout,*) "Error reading electrostatic energy parameters."
1050       goto 999
1051   117 write (iout,*) "Error reading side chain interaction parameters."
1052       goto 999
1053   118 write (iout,*) "Error reading SCp interaction parameters."
1054       goto 999
1055   119 write (iout,*) "Error reading SCCOR parameters"
1056   999 continue
1057 #ifdef MPI
1058       call MPI_Finalize(Ierror)
1059 #endif
1060       stop
1061       return
1062       end
1063
1064
1065       subroutine getenv_loc(var, val)
1066       character(*) var, val
1067
1068 #ifdef WINIFL
1069       character(2000) line
1070       external ilen
1071
1072       open (196,file='env',status='old',readonly,shared)
1073       iread=0
1074 c      write(*,*)'looking for ',var
1075 10    read(196,*,err=11,end=11)line
1076       iread=index(line,var)
1077 c      write(*,*)iread,' ',var,' ',line
1078       if (iread.eq.0) go to 10 
1079 c      write(*,*)'---> ',line
1080 11    continue
1081       if(iread.eq.0) then
1082 c       write(*,*)'CHUJ'
1083        val=''
1084       else
1085        iread=iread+ilen(var)+1
1086        read (line(iread:),*,err=12,end=12) val
1087 c       write(*,*)'OK: ',var,' = ',val
1088       endif
1089       close(196)
1090       return
1091 12    val=''
1092       close(196)
1093 #elif (defined CRAY)
1094       integer lennam,lenval,ierror
1095 c
1096 c        getenv using a POSIX call, useful on the T3D
1097 c        Sept 1996, comment out error check on advice of H. Pritchard
1098 c
1099       lennam = len(var)
1100       if(lennam.le.0) stop '--error calling getenv--'
1101       call pxfgetenv(var,lennam,val,lenval,ierror)
1102 c-HP- if(ierror.ne.0) stop '--error returned by pxfgetenv--'
1103 #else
1104       call getenv(var,val)
1105 #endif
1106
1107       return
1108       end