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