Merge branch 'master' of mmka:unres into devel
[unres.git] / source / wham / src / parmread.F
1       subroutine parmread(iparm,*)
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       implicit real*8 (a-h,o-z)
7       include 'DIMENSIONS'
8       include 'DIMENSIONS.ZSCOPT'
9       include 'DIMENSIONS.FREE'
10       include 'COMMON.IOUNITS'
11       include 'COMMON.CHAIN'
12       include 'COMMON.INTERACT'
13       include 'COMMON.GEO'
14       include 'COMMON.LOCAL'
15       include 'COMMON.TORSION'
16       include 'COMMON.FFIELD'
17       include 'COMMON.NAMES'
18       include 'COMMON.SBRIDGE'
19       include 'COMMON.WEIGHTS'
20       include 'COMMON.ENEPS'
21       include 'COMMON.SCCOR'
22       include 'COMMON.SCROT'
23       include 'COMMON.FREE'
24       character*1 t1,t2,t3
25       character*1 onelett(4) /"G","A","P","D"/
26       logical lprint
27       dimension blower(3,3,maxlob)
28       character*800 controlcard
29       character*256 bondname_t,thetname_t,rotname_t,torname_t,
30      & tordname_t,fouriername_t,elename_t,sidename_t,scpname_t,
31      & sccorname_t
32       integer ilen
33       external ilen
34       character*16 key
35       integer iparm
36       double precision ip,mp
37 C
38 C Body
39 C
40 C Set LPRINT=.TRUE. for debugging
41       dwa16=2.0d0**(1.0d0/6.0d0)
42       lprint=.false.
43       itypro=20
44 C Assign virtual-bond length
45       vbl=3.8D0
46       vblinv=1.0D0/vbl
47       vblinv2=vblinv*vblinv
48       call card_concat(controlcard,.true.)
49       wname(4)="WCORRH"
50       do i=1,n_ene
51         key = wname(i)(:ilen(wname(i)))
52         call reada(controlcard,key(:ilen(key)),ww(i),1.0d0)
53       enddo
54
55       write (iout,*) "iparm",iparm," myparm",myparm
56 c If reading not own parameters, skip assignment
57
58       if (iparm.eq.myparm .or. .not.separate_parset) then
59
60 c
61 c Setup weights for UNRES
62 c
63       wsc=ww(1)
64       wscp=ww(2)
65       welec=ww(3)
66       wcorr=ww(4)
67       wcorr5=ww(5)
68       wcorr6=ww(6)
69       wel_loc=ww(7)
70       wturn3=ww(8)
71       wturn4=ww(9)
72       wturn6=ww(10)
73       wang=ww(11)
74       wscloc=ww(12)
75       wtor=ww(13)
76       wtor_d=ww(14)
77       wvdwpp=ww(16)
78       wstrain=ww(15)
79       wbond=ww(18)
80       wsccor=ww(19)
81
82       endif
83
84       call card_concat(controlcard,.false.)
85
86 c Return if not own parameters
87
88       if (iparm.ne.myparm .and. separate_parset) return
89
90       call reads(controlcard,"BONDPAR",bondname_t,bondname)
91       open (ibond,file=bondname_t,status='old')
92       rewind(ibond)
93       call reads(controlcard,"THETPAR",thetname_t,thetname)
94       open (ithep,file=thetname_t,status='old')
95       rewind(ithep) 
96       call reads(controlcard,"ROTPAR",rotname_t,rotname)
97       open (irotam,file=rotname_t,status='old')
98       rewind(irotam)
99       call reads(controlcard,"TORPAR",torname_t,torname)
100       open (itorp,file=torname_t,status='old')
101       rewind(itorp)
102       call reads(controlcard,"TORDPAR",tordname_t,tordname)
103       open (itordp,file=tordname_t,status='old')
104       rewind(itordp)
105       call reads(controlcard,"SCCORAR",sccorname_t,sccorname)
106       open (isccor,file=sccorname_t,status='old')
107       rewind(isccor)
108       call reads(controlcard,"FOURIER",fouriername_t,fouriername)
109       open (ifourier,file=fouriername_t,status='old')
110       rewind(ifourier)
111       call reads(controlcard,"ELEPAR",elename_t,elename)
112       open (ielep,file=elename_t,status='old')
113       rewind(ielep)
114       call reads(controlcard,"SIDEPAR",sidename_t,sidename)
115       open (isidep,file=sidename_t,status='old')
116       rewind(isidep)
117       call reads(controlcard,"SCPPAR",scpname_t,scpname)
118       open (iscpp,file=scpname_t,status='old')
119       rewind(iscpp)
120       write (iout,*) "Parameter set:",iparm
121       write (iout,*) "Energy-term weights:"
122       do i=1,n_ene
123         write (iout,'(a16,f10.5)') wname(i),ww(i)
124       enddo
125       write (iout,*) "Sidechain potential file        : ",
126      &  sidename_t(:ilen(sidename_t))
127 #ifndef OLDSCP
128       write (iout,*) "SCp potential file              : ",
129      &  scpname_t(:ilen(scpname_t))
130 #endif  
131       write (iout,*) "Electrostatic potential file    : ",
132      &  elename_t(:ilen(elename_t))
133       write (iout,*) "Cumulant coefficient file       : ",
134      &  fouriername_t(:ilen(fouriername_t))
135       write (iout,*) "Torsional parameter file        : ",
136      &  torname_t(:ilen(torname_t))
137       write (iout,*) "Double torsional parameter file : ",
138      &  tordname_t(:ilen(tordname_t))
139       write (iout,*) "Backbone-rotamer parameter file : ",
140      &  sccorname(:ilen(sccorname))
141       write (iout,*) "Bond & inertia constant file    : ",
142      &  bondname_t(:ilen(bondname_t))
143       write (iout,*) "Bending parameter file          : ",
144      &  thetname_t(:ilen(thetname_t))
145       write (iout,*) "Rotamer parameter file          : ",
146      &  rotname_t(:ilen(rotname_t))
147
148 c
149 c Read the virtual-bond parameters, masses, and moments of inertia
150 c and Stokes' radii of the peptide group and side chains
151 c
152 #ifdef CRYST_BOND
153       read (ibond,*) vbldp0,akp
154       do i=1,ntyp
155         nbondterm(i)=1
156         read (ibond,*) vbldsc0(1,i),aksc(1,i)
157         dsc(i) = vbldsc0(1,i)
158         if (i.eq.10) then
159           dsc_inv(i)=0.0D0
160         else
161           dsc_inv(i)=1.0D0/dsc(i)
162         endif
163       enddo
164 #else
165       read (ibond,*) ijunk,vbldp0,akp,rjunk
166       do i=1,ntyp
167         read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),
168      &   j=1,nbondterm(i))
169         dsc(i) = vbldsc0(1,i)
170         if (i.eq.10) then
171           dsc_inv(i)=0.0D0
172         else
173           dsc_inv(i)=1.0D0/dsc(i)
174         endif
175       enddo
176 #endif
177       if (lprint) then
178         write(iout,'(/a/)')"Force constants virtual bonds:"
179         write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K',
180      &   'inertia','Pstok'
181         write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0
182         do i=1,ntyp
183           write (iout,'(a10,i3,6f10.5)') restyp(i),nbondterm(i),
184      &      vbldsc0(1,i),aksc(1,i),abond0(1,i)
185           do j=2,nbondterm(i)
186             write (iout,'(13x,3f10.5)')
187      &        vbldsc0(j,i),aksc(j,i),abond0(j,i)
188           enddo
189         enddo
190       endif
191 #ifdef CRYST_THETA
192 C
193 C Read the parameters of the probability distribution/energy expression 
194 C of the virtual-bond valence angles theta
195 C
196       do i=1,ntyp
197         read (ithep,*) a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
198         read (ithep,*) (polthet(j,i),j=0,3)
199         read (ithep,*) (gthet(j,i),j=1,3)
200         read (ithep,*) theta0(i),sig0(i),sigc0(i)
201         sigc0(i)=sigc0(i)**2
202       enddo
203       close (ithep)
204       if (lprint) then
205 c       write (iout,'(a)') 
206 c    &    'Parameters of the virtual-bond valence angles:'
207 c       write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',
208 c    & '    ATHETA0   ','         A1   ','        A2    ',
209 c    & '        B1    ','         B2   '        
210 c       do i=1,ntyp
211 c         write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
212 c    &        a0thet(i),(athet(j,i),j=1,2),(bthet(j,i),j=1,2)
213 c       enddo
214 c       write (iout,'(/a/9x,5a/79(1h-))') 
215 c    & 'Parameters of the expression for sigma(theta_c):',
216 c    & '     ALPH0    ','      ALPH1   ','     ALPH2    ',
217 c    & '     ALPH3    ','    SIGMA0C   '        
218 c       do i=1,ntyp
219 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,
220 c    &      (polthet(j,i),j=0,3),sigc0(i) 
221 c       enddo
222 c       write (iout,'(/a/9x,5a/79(1h-))') 
223 c    & 'Parameters of the second gaussian:',
224 c    & '    THETA0    ','     SIGMA0   ','        G1    ',
225 c    & '        G2    ','         G3   '        
226 c       do i=1,ntyp
227 c         write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i),i,theta0(i),
228 c    &       sig0(i),(gthet(j,i),j=1,3)
229 c       enddo
230         write (iout,'(a)') 
231      &    'Parameters of the virtual-bond valence angles:'
232         write (iout,'(/a/9x,5a/79(1h-))') 
233      & 'Coefficients of expansion',
234      & '     theta0   ','    a1*10^2   ','   a2*10^2    ',
235      & '   b1*10^1    ','    b2*10^1   '        
236         do i=1,ntyp
237           write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),
238      &        a0thet(i),(100*athet(j,i),j=1,2),(10*bthet(j,i),j=1,2)
239         enddo
240         write (iout,'(/a/9x,5a/79(1h-))') 
241      & 'Parameters of the expression for sigma(theta_c):',
242      & ' alpha0       ','  alph1       ',' alph2        ',
243      & ' alhp3        ','   sigma0c    '        
244         do i=1,ntyp
245           write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i),
246      &      (polthet(j,i),j=0,3),sigc0(i) 
247         enddo
248         write (iout,'(/a/9x,5a/79(1h-))') 
249      & 'Parameters of the second gaussian:',
250      & '    theta0    ','  sigma0*10^2 ','      G1*10^-1',
251      & '        G2    ','   G3*10^1    '        
252         do i=1,ntyp
253           write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i),theta0(i),
254      &       100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
255         enddo
256       endif
257 #else
258 C
259 C Read the parameters of Utheta determined from ab initio surfaces
260 C Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
261 C
262       read (ithep,*) nthetyp,ntheterm,ntheterm2,
263      &  ntheterm3,nsingle,ndouble
264       nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
265       read (ithep,*) (ithetyp(i),i=1,ntyp1)
266       do i=1,maxthetyp
267         do j=1,maxthetyp
268           do k=1,maxthetyp
269             aa0thet(i,j,k)=0.0d0
270             do l=1,ntheterm
271               aathet(l,i,j,k)=0.0d0
272             enddo
273             do l=1,ntheterm2
274               do m=1,nsingle
275                 bbthet(m,l,i,j,k)=0.0d0
276                 ccthet(m,l,i,j,k)=0.0d0
277                 ddthet(m,l,i,j,k)=0.0d0
278                 eethet(m,l,i,j,k)=0.0d0
279               enddo
280             enddo
281             do l=1,ntheterm3
282               do m=1,ndouble
283                 do mm=1,ndouble
284                  ffthet(mm,m,l,i,j,k)=0.0d0
285                  ggthet(mm,m,l,i,j,k)=0.0d0
286                 enddo
287               enddo
288             enddo
289           enddo
290         enddo
291       enddo
292       do i=1,nthetyp
293         do j=1,nthetyp
294           do k=1,nthetyp
295             read (ithep,'(3a)') res1,res2,res3
296             read (ithep,*) aa0thet(i,j,k)
297             read (ithep,*)(aathet(l,i,j,k),l=1,ntheterm)
298             read (ithep,*)
299      &       ((bbthet(lll,ll,i,j,k),lll=1,nsingle),
300      &        (ccthet(lll,ll,i,j,k),lll=1,nsingle),
301      &        (ddthet(lll,ll,i,j,k),lll=1,nsingle),
302      &        (eethet(lll,ll,i,j,k),lll=1,nsingle),ll=1,ntheterm2)
303             read (ithep,*)
304      &      (((ffthet(llll,lll,ll,i,j,k),ffthet(lll,llll,ll,i,j,k),
305      &         ggthet(llll,lll,ll,i,j,k),ggthet(lll,llll,ll,i,j,k),
306      &         llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
307           enddo
308         enddo
309       enddo
310 C
311 C For dummy ends assign glycine-type coefficients of theta-only terms; the
312 C coefficients of theta-and-gamma-dependent terms are zero.
313 C
314       do i=1,nthetyp
315         do j=1,nthetyp
316           do l=1,ntheterm
317             aathet(l,i,j,nthetyp+1)=aathet(l,i,j,1)
318             aathet(l,nthetyp+1,i,j)=aathet(l,1,i,j)
319           enddo
320           aa0thet(i,j,nthetyp+1)=aa0thet(i,j,1)
321           aa0thet(nthetyp+1,i,j)=aa0thet(1,i,j)
322         enddo
323         do l=1,ntheterm
324           aathet(l,nthetyp+1,i,nthetyp+1)=aathet(l,1,i,1)
325         enddo
326         aa0thet(nthetyp+1,i,nthetyp+1)=aa0thet(1,i,1)
327       enddo
328 C
329 C Control printout of the coefficients of virtual-bond-angle potentials
330 C
331       if (lprint) then
332         write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
333         do i=1,nthetyp+1
334           do j=1,nthetyp+1
335             do k=1,nthetyp+1
336               write (iout,'(//4a)')
337      &         'Type ',onelett(i),onelett(j),onelett(k)
338               write (iout,'(//a,10x,a)') " l","a[l]"
339               write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k)
340               write (iout,'(i2,1pe15.5)')
341      &           (l,aathet(l,i,j,k),l=1,ntheterm)
342             do l=1,ntheterm2
343               write (iout,'(//2h m,4(9x,a,3h[m,i1,1h]))')
344      &          "b",l,"c",l,"d",l,"e",l
345               do m=1,nsingle
346                 write (iout,'(i2,4(1pe15.5))') m,
347      &          bbthet(m,l,i,j,k),ccthet(m,l,i,j,k),
348      &          ddthet(m,l,i,j,k),eethet(m,l,i,j,k)
349               enddo
350             enddo
351             do l=1,ntheterm3
352               write (iout,'(//3hm,n,4(6x,a,5h[m,n,i1,1h]))')
353      &          "f+",l,"f-",l,"g+",l,"g-",l
354               do m=2,ndouble
355                 do n=1,m-1
356                   write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,
357      &              ffthet(n,m,l,i,j,k),ffthet(m,n,l,i,j,k),
358      &              ggthet(n,m,l,i,j,k),ggthet(m,n,l,i,j,k)
359                 enddo
360               enddo
361             enddo
362           enddo
363         enddo
364       enddo
365       call flush(iout)
366       endif
367 #endif
368
369 #ifdef CRYST_SC
370 C
371 C Read the parameters of the probability distribution/energy expression
372 C of the side chains.
373 C
374       do i=1,ntyp
375         read (irotam,'(3x,i3,f8.3)') nlob(i),dsc(i)
376         if (i.eq.10) then
377           dsc_inv(i)=0.0D0
378         else
379           dsc_inv(i)=1.0D0/dsc(i)
380         endif
381         if (i.ne.10) then
382         do j=1,nlob(i)
383           do k=1,3
384             do l=1,3
385               blower(l,k,j)=0.0D0
386             enddo
387           enddo
388         enddo  
389         bsc(1,i)=0.0D0
390         read(irotam,*)(censc(k,1,i),k=1,3),((blower(k,l,1),l=1,k),k=1,3)
391         do j=2,nlob(i)
392           read (irotam,*) bsc(j,i)
393           read (irotam,*) (censc(k,j,i),k=1,3),
394      &                                 ((blower(k,l,j),l=1,k),k=1,3)
395         enddo
396         do j=1,nlob(i)
397           do k=1,3
398             do l=1,k
399               akl=0.0D0
400               do m=1,3
401                 akl=akl+blower(k,m,j)*blower(l,m,j)
402               enddo
403               gaussc(k,l,j,i)=akl
404               gaussc(l,k,j,i)=akl
405             enddo
406           enddo 
407         enddo
408         endif
409       enddo
410       close (irotam)
411       if (lprint) then
412         write (iout,'(/a)') 'Parameters of side-chain local geometry'
413         do i=1,ntyp
414           nlobi=nlob(i)
415           if (nlobi.gt.0) then
416           write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i),
417      &     ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
418 c          write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
419 c          write (iout,'(a,f10.4,4(16x,f10.4))')
420 c     &                             'Center  ',(bsc(j,i),j=1,nlobi)
421 c          write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),j=1,nlobi)
422            write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))')
423      &                             'log h',(bsc(j,i),j=1,nlobi)
424            write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))')
425      &    'x',((censc(k,j,i),k=1,3),j=1,nlobi)
426 c          write (iout,'(a)')
427 c         do j=1,nlobi
428 c           ind=0
429 c           do k=1,3
430 c             do l=1,k
431 c              ind=ind+1
432 c              blower(k,l,j)=gaussc(ind,j,i)
433 c             enddo
434 c           enddo
435 c         enddo
436           do k=1,3
437             write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))')
438      &                 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
439           enddo
440           endif
441         enddo
442       endif
443 #else
444 C
445 C Read scrot parameters for potentials determined from all-atom AM1 calculations
446 C added by Urszula Kozlowska 07/11/2007
447 C
448       do i=1,ntyp
449         read (irotam,*)
450        if (i.eq.10) then
451          read (irotam,*)
452        else
453          do j=1,65
454            read(irotam,*) sc_parmin(j,i)
455          enddo
456        endif
457       enddo
458 #endif
459       close(irotam)
460 #ifdef CRYST_TOR
461 C
462 C Read torsional parameters in old format
463 C
464       read (itorp,*) ntortyp,nterm_old
465       write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
466       read (itorp,*) (itortyp(i),i=1,ntyp)
467       do i=1,ntortyp
468         do j=1,ntortyp
469           read (itorp,'(a)')
470           do k=1,nterm_old
471             read (itorp,*) kk,v1(k,j,i),v2(k,j,i) 
472           enddo
473         enddo
474       enddo
475       close (itorp)
476       if (lprint) then
477         write (iout,'(/a/)') 'Torsional constants:'
478         do i=1,ntortyp
479           do j=1,ntortyp
480             write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
481             write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
482           enddo
483         enddo
484       endif
485
486
487 #else
488 C
489 C Read torsional parameters
490 C
491       read (itorp,*) ntortyp
492       read (itorp,*) (itortyp(i),i=1,ntyp)
493       write (iout,*) 'ntortyp',ntortyp
494       do i=1,ntortyp
495         do j=1,ntortyp
496           read (itorp,*) nterm(i,j),nlor(i,j)
497           v0ij=0.0d0
498           si=-1.0d0
499           do k=1,nterm(i,j)
500             read (itorp,*) kk,v1(k,i,j),v2(k,i,j) 
501             v0ij=v0ij+si*v1(k,i,j)
502             si=-si
503           enddo
504           do k=1,nlor(i,j)
505             read (itorp,*) kk,vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j) 
506             v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
507           enddo
508           v0(i,j)=v0ij
509         enddo
510       enddo
511       close (itorp)
512       if (lprint) then
513         write (iout,'(/a/)') 'Torsional constants:'
514         do i=1,ntortyp
515           do j=1,ntortyp
516             write (iout,*) 'ityp',i,' jtyp',j
517             write (iout,*) 'Fourier constants'
518             do k=1,nterm(i,j)
519               write (iout,'(2(1pe15.5))') v1(k,i,j),v2(k,i,j)
520             enddo
521             write (iout,*) 'Lorenz constants'
522             do k=1,nlor(i,j)
523               write (iout,'(3(1pe15.5))') 
524      &         vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
525             enddo
526           enddo
527         enddo
528       endif
529 C
530 C 6/23/01 Read parameters for double torsionals
531 C
532       do i=1,ntortyp
533         do j=1,ntortyp
534           do k=1,ntortyp
535             read (itordp,'(3a1)') t1,t2,t3
536             if (t1.ne.onelett(i) .or. t2.ne.onelett(j) 
537      &        .or. t3.ne.onelett(k)) then
538               write (iout,*) "Error in double torsional parameter file",
539      &         i,j,k,t1,t2,t3
540                stop "Error in double torsional parameter file"
541             endif
542             read (itordp,*) ntermd_1(i,j,k),ntermd_2(i,j,k)
543             read (itordp,*) (v1c(1,l,i,j,k),l=1,ntermd_1(i,j,k))
544             read (itordp,*) (v1s(1,l,i,j,k),l=1,ntermd_1(i,j,k))
545             read (itordp,*) (v1c(2,l,i,j,k),l=1,ntermd_1(i,j,k))
546             read (itordp,*) (v1s(2,l,i,j,k),l=1,ntermd_1(i,j,k))
547             read (itordp,*) ((v2c(l,m,i,j,k),v2c(m,l,i,j,k),
548      &       v2s(l,m,i,j,k),v2s(m,l,i,j,k),m=1,l-1),l=1,ntermd_2(i,j,k))
549           enddo
550         enddo
551       enddo
552       if (lprint) then
553       write (iout,*) 
554       write (iout,*) 'Constants for double torsionals'
555       do i=1,ntortyp
556         do j=1,ntortyp 
557           do k=1,ntortyp
558             write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,
559      &        ' nsingle',ntermd_1(i,j,k),' ndouble',ntermd_2(i,j,k)
560             write (iout,*)
561             write (iout,*) 'Single angles:'
562             do l=1,ntermd_1(i,j,k)
563               write (iout,'(i5,2f10.5,5x,2f10.5)') l,
564      &           v1c(1,l,i,j,k),v1s(1,l,i,j,k),
565      &           v1c(2,l,i,j,k),v1s(2,l,i,j,k)
566             enddo
567             write (iout,*)
568             write (iout,*) 'Pairs of angles:'
569             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
570             do l=1,ntermd_2(i,j,k)
571               write (iout,'(i5,20f10.5)') 
572      &         l,(v2c(l,m,i,j,k),m=1,ntermd_2(i,j,k))
573             enddo
574             write (iout,*)
575             write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k))
576             do l=1,ntermd_2(i,j,k)
577               write (iout,'(i5,20f10.5)') 
578      &         l,(v2s(l,m,i,j,k),m=1,ntermd_2(i,j,k))
579             enddo
580             write (iout,*)
581           enddo
582         enddo
583       enddo
584       endif
585 #endif
586 C Read of Side-chain backbone correlation parameters
587 C Modified 11 May 2012 by Adasko
588 CCC
589 C
590       read (isccor,*) nsccortyp
591       read (isccor,*) (isccortyp(i),i=1,ntyp)
592 c      write (iout,*) 'ntortyp',ntortyp
593       maxinter=3
594 cc maxinter is maximum interaction sites
595       do l=1,maxinter
596       do i=1,nsccortyp
597         do j=1,nsccortyp
598           read (isccor,*) nterm_sccor(i,j),nlor_sccor(i,j)
599           v0ijsccor=0.0d0
600           si=-1.0d0
601
602           do k=1,nterm_sccor(i,j)
603             read (isccor,*) kk,v1sccor(k,l,i,j)
604      &    ,v2sccor(k,l,i,j)
605             v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
606             si=-si
607           enddo
608           do k=1,nlor_sccor(i,j)
609             read (isccor,*) kk,vlor1sccor(k,i,j),
610      &        vlor2sccor(k,i,j),vlor3sccor(k,i,j)
611             v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/
612      &(1+vlor3sccor(k,i,j)**2)
613           enddo
614           v0sccor(i,j)=v0ijsccor
615         enddo
616       enddo
617       enddo
618       close (isccor)
619
620       if (lprint) then
621         write (iout,'(/a/)') 'Torsional constants:'
622         do i=1,nsccortyp
623           do j=1,nsccortyp
624             write (iout,*) 'ityp',i,' jtyp',j
625             write (iout,*) 'Fourier constants'
626             do k=1,nterm_sccor(i,j)
627            write (iout,'(2(1pe15.5))')v1sccor(k,l,i,j),v2sccor(k,l,i,j)
628             enddo
629             write (iout,*) 'Lorenz constants'
630             do k=1,nlor_sccor(i,j)
631               write (iout,'(3(1pe15.5))')
632      &         vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
633             enddo
634           enddo
635         enddo
636       endif
637
638 C
639 C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
640 C         interaction energy of the Gly, Ala, and Pro prototypes.
641 C
642       read (ifourier,*) nloctyp
643       do i=1,nloctyp
644         read (ifourier,*)
645         read (ifourier,*) (b(ii,i),ii=1,13)
646         if (lprint) then
647         write (iout,*) 'Type',i
648         write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
649         endif
650         B1(1,i)  = b(3,i)
651         B1(2,i)  = b(5,i)
652         B1tilde(1,i) = b(3,i)
653         B1tilde(2,i) =-b(5,i) 
654         B2(1,i)  = b(2,i)
655         B2(2,i)  = b(4,i)
656         CC(1,1,i)= b(7,i)
657         CC(2,2,i)=-b(7,i)
658         CC(2,1,i)= b(9,i)
659         CC(1,2,i)= b(9,i)
660         Ctilde(1,1,i)=b(7,i)
661         Ctilde(1,2,i)=b(9,i)
662         Ctilde(2,1,i)=-b(9,i)
663         Ctilde(2,2,i)=b(7,i)
664         DD(1,1,i)= b(6,i)
665         DD(2,2,i)=-b(6,i)
666         DD(2,1,i)= b(8,i)
667         DD(1,2,i)= b(8,i)
668         Dtilde(1,1,i)=b(6,i)
669         Dtilde(1,2,i)=b(8,i)
670         Dtilde(2,1,i)=-b(8,i)
671         Dtilde(2,2,i)=b(6,i)
672         EE(1,1,i)= b(10,i)+b(11,i)
673         EE(2,2,i)=-b(10,i)+b(11,i)
674         EE(2,1,i)= b(12,i)-b(13,i)
675         EE(1,2,i)= b(12,i)+b(13,i)
676       enddo
677       if (lprint) then
678       do i=1,nloctyp
679         write (iout,*) 'Type',i
680         write (iout,*) 'B1'
681 c        write (iout,'(f10.5)') B1(:,i)
682         write(iout,*) B1(1,i),B1(2,i)
683         write (iout,*) 'B2'
684 c        write (iout,'(f10.5)') B2(:,i)
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,'(/a)') 'Electrostatic interaction constants:'
705         write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') 
706      &            'IT','JT','APP','BPP','AEL6','AEL3'
707       endif
708       read (ielep,*) ((epp(i,j),j=1,2),i=1,2)
709       read (ielep,*) ((rpp(i,j),j=1,2),i=1,2)
710       read (ielep,*) ((elpp6(i,j),j=1,2),i=1,2)
711       read (ielep,*) ((elpp3(i,j),j=1,2),i=1,2)
712       close (ielep)
713       do i=1,2
714         do j=1,2
715         rri=rpp(i,j)**6
716         app (i,j)=epp(i,j)*rri*rri 
717         bpp (i,j)=-2.0D0*epp(i,j)*rri
718         ael6(i,j)=elpp6(i,j)*4.2D0**6
719         ael3(i,j)=elpp3(i,j)*4.2D0**3
720         if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),
721      &                    ael6(i,j),ael3(i,j)
722         enddo
723       enddo
724 C
725 C Read side-chain interaction parameters.
726 C
727       read (isidep,*) ipot,expon
728       if (ipot.lt.1 .or. ipot.gt.5) then
729         write (iout,'(2a)') 'Error while reading SC interaction',
730      &               'potential file - unknown potential type.'
731         stop
732       endif
733       expon2=expon/2
734       write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),
735      & ', exponents are ',expon,2*expon 
736       goto (10,20,30,30,40) ipot
737 C----------------------- LJ potential ---------------------------------
738    10 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),(sigma0(i),i=1,ntyp)
739       if (lprint) then
740         write (iout,'(/a/)') 'Parameters of the LJ potential:'
741         write (iout,'(a/)') 'The epsilon array:'
742         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
743         write (iout,'(/a)') 'One-body parameters:'
744         write (iout,'(a,4x,a)') 'residue','sigma'
745         write (iout,'(a3,6x,f10.5)') (restyp(i),sigma0(i),i=1,ntyp)
746       endif
747       goto 50
748 C----------------------- LJK potential --------------------------------
749    20 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
750      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
751       if (lprint) then
752         write (iout,'(/a/)') 'Parameters of the LJK potential:'
753         write (iout,'(a/)') 'The epsilon array:'
754         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
755         write (iout,'(/a)') 'One-body parameters:'
756         write (iout,'(a,4x,2a)') 'residue','   sigma  ','    r0    '
757         write (iout,'(a3,6x,2f10.5)') (restyp(i),sigma0(i),rr0(i),
758      &        i=1,ntyp)
759       endif
760       goto 50
761 C---------------------- GB or BP potential -----------------------------
762    30 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
763      &  (sigma0(i),i=1,ntyp),(sigii(i),i=1,ntyp),(chip0(i),i=1,ntyp),
764      &  (alp(i),i=1,ntyp)
765 C For the GB potential convert sigma'**2 into chi'
766       if (ipot.eq.4) then
767         do i=1,ntyp
768           chip(i)=(chip0(i)-1.0D0)/(chip0(i)+1.0D0)
769         enddo
770       endif
771       if (lprint) then
772         write (iout,'(/a/)') 'Parameters of the BP potential:'
773         write (iout,'(a/)') 'The epsilon array:'
774         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
775         write (iout,'(/a)') 'One-body parameters:'
776         write (iout,'(a,4x,4a)') 'residue','   sigma  ','s||/s_|_^2',
777      &       '    chip  ','    alph  '
778         write (iout,'(a3,6x,4f10.5)') (restyp(i),sigma0(i),sigii(i),
779      &                     chip(i),alp(i),i=1,ntyp)
780       endif
781       goto 50
782 C--------------------- GBV potential -----------------------------------
783    40 read (isidep,*)((eps(i,j),j=i,ntyp),i=1,ntyp),
784      &  (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),
785      &  (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
786       if (lprint) then
787         write (iout,'(/a/)') 'Parameters of the GBV potential:'
788         write (iout,'(a/)') 'The epsilon array:'
789         call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
790         write (iout,'(/a)') 'One-body parameters:'
791         write (iout,'(a,4x,5a)') 'residue','   sigma  ','    r0    ',
792      &      's||/s_|_^2','    chip  ','    alph  '
793         write (iout,'(a3,6x,5f10.5)') (restyp(i),sigma0(i),rr0(i),
794      &           sigii(i),chip(i),alp(i),i=1,ntyp)
795       endif
796    50 continue
797       close (isidep)
798 C-----------------------------------------------------------------------
799 C Calculate the "working" parameters of SC interactions.
800       do i=2,ntyp
801         do j=1,i-1
802           eps(i,j)=eps(j,i)
803         enddo
804       enddo
805       do i=1,ntyp
806         do j=i,ntyp
807           sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
808           sigma(j,i)=sigma(i,j)
809           rs0(i,j)=dwa16*sigma(i,j)
810           rs0(j,i)=rs0(i,j)
811         enddo
812       enddo
813       if (lprint) write (iout,'(/a/10x,7a/72(1h-))') 
814      & 'Working parameters of the SC interactions:',
815      & '     a    ','     b    ','   augm   ','  sigma ','   r0   ',
816      & '  chi1   ','   chi2   ' 
817       do i=1,ntyp
818         do j=i,ntyp
819           epsij=eps(i,j)
820           if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
821             rrij=sigma(i,j)
822           else
823             rrij=rr0(i)+rr0(j)
824           endif
825           r0(i,j)=rrij
826           r0(j,i)=rrij
827           rrij=rrij**expon
828           epsij=eps(i,j)
829           sigeps=dsign(1.0D0,epsij)
830           epsij=dabs(epsij)
831           aa(i,j)=epsij*rrij*rrij
832           bb(i,j)=-sigeps*epsij*rrij
833           aa(j,i)=aa(i,j)
834           bb(j,i)=bb(i,j)
835           if (ipot.gt.2) then
836             sigt1sq=sigma0(i)**2
837             sigt2sq=sigma0(j)**2
838             sigii1=sigii(i)
839             sigii2=sigii(j)
840             ratsig1=sigt2sq/sigt1sq
841             ratsig2=1.0D0/ratsig1
842             chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
843             if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
844             rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
845           else
846             rsum_max=sigma(i,j)
847           endif
848 c         if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
849             sigmaii(i,j)=rsum_max
850             sigmaii(j,i)=rsum_max 
851 c         else
852 c           sigmaii(i,j)=r0(i,j)
853 c           sigmaii(j,i)=r0(i,j)
854 c         endif
855 cd        write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
856           if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
857             r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
858             augm(i,j)=epsij*r_augm**(2*expon)
859 c           augm(i,j)=0.5D0**(2*expon)*aa(i,j)
860             augm(j,i)=augm(i,j)
861           else
862             augm(i,j)=0.0D0
863             augm(j,i)=0.0D0
864           endif
865           if (lprint) then
866             write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') 
867      &      restyp(i),restyp(j),aa(i,j),bb(i,j),augm(i,j),
868      &      sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
869           endif
870         enddo
871       enddo
872 C
873 C Define the SC-p interaction constants
874 C
875 #ifdef OLDSCP
876       do i=1,20
877 C "Soft" SC-p repulsion (causes helices to be too flat, but facilitates 
878 C helix formation)
879 c       aad(i,1)=0.3D0*4.0D0**12
880 C Following line for constants currently implemented
881 C "Hard" SC-p repulsion (gives correct turn spacing in helices)
882         aad(i,1)=1.5D0*4.0D0**12
883 c       aad(i,1)=0.17D0*5.6D0**12
884         aad(i,2)=aad(i,1)
885 C "Soft" SC-p repulsion
886         bad(i,1)=0.0D0
887 C Following line for constants currently implemented
888 c       aad(i,1)=0.3D0*4.0D0**6
889 C "Hard" SC-p repulsion
890         bad(i,1)=3.0D0*4.0D0**6
891 c       bad(i,1)=-2.0D0*0.17D0*5.6D0**6
892         bad(i,2)=bad(i,1)
893 c       aad(i,1)=0.0D0
894 c       aad(i,2)=0.0D0
895 c       bad(i,1)=1228.8D0
896 c       bad(i,2)=1228.8D0
897       enddo
898 #else
899 C
900 C 8/9/01 Read the SC-p interaction constants from file
901 C
902       do i=1,ntyp
903         read (iscpp,*) (eps_scp(i,j),rscp(i,j),j=1,2)
904       enddo
905       do i=1,ntyp
906         aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
907         aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
908         bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
909         bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
910       enddo
911
912       if (lprint) then
913         write (iout,*) "Parameters of SC-p interactions:"
914         do i=1,20
915           write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),
916      &     eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
917         enddo
918       endif
919 #endif
920 C
921 C Define the constants of the disulfide bridge
922 C
923       ebr=-5.50D0
924 c
925 c Old arbitrary potential - commented out.
926 c
927 c      dbr= 4.20D0
928 c      fbr= 3.30D0
929 c
930 c Constants of the disulfide-bond potential determined based on the RHF/6-31G**
931 c energy surface of diethyl disulfide.
932 c A. Liwo and U. Kozlowska, 11/24/03
933 c
934       D0CM = 3.78d0
935       AKCM = 15.1d0
936       AKTH = 11.0d0
937       AKCT = 12.0d0
938       V1SS =-1.08d0
939       V2SS = 7.61d0
940       V3SS = 13.7d0
941
942       if (lprint) then
943       write (iout,'(/a)') "Disulfide bridge parameters:"
944       write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
945       write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
946       write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
947       write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,
948      & ' v3ss:',v3ss
949       endif
950       return
951       end