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