8 use control_data, only:maxterm_sccor
10 !-----------------------------------------------------------------------------
11 ! Max. number of residue types and parameters in expressions for
12 ! virtual-bond angle bending potentials
13 ! integer,parameter :: maxthetyp=3
14 ! integer,parameter :: maxthetyp1=maxthetyp+1
16 ! & maxtheterm2=6,maxtheterm3=4,maxsingle=6,maxdouble=4,
17 ! & mmaxtheterm=maxtheterm)
18 !-----------------------------------------------------------------------------
19 ! Max. number of types of dihedral angles & multiplicity of torsional barriers
20 ! and the number of terms in double torsionals
21 ! integer,parameter :: maxlor=3,maxtermd_1=8,maxtermd_2=8
22 ! parameter (maxtor=4,maxterm=10)
23 !-----------------------------------------------------------------------------
24 ! Max number of torsional terms in SCCOR
25 !el integer,parameter :: maxterm_sccor=6
26 !-----------------------------------------------------------------------------
27 character(len=1),dimension(:),allocatable :: secstruc !(maxres)
28 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
33 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
34 !-----------------------------------------------------------------------------
36 !-----------------------------------------------------------------------------
37 subroutine write_rbank(jlee,adif,nft)
40 use geometry_data, only: nres,rad2deg
41 ! implicit real*8 (a-h,o-z)
42 ! include 'DIMENSIONS'
43 ! include 'COMMON.IOUNITS'
44 ! include 'COMMON.CSA'
45 ! include 'COMMON.BANK'
46 ! include 'COMMON.CHAIN'
47 ! include 'COMMON.GEO'
49 integer :: nft,i,k,j,l,jlee
52 open(icsa_rbank,file=csa_rbank,status="unknown")
53 write (icsa_rbank,900) jlee,nbank,nstep,nft,icycle,adif
55 write (icsa_rbank,952) k,rene(k),rrmsn(k),rpncn(k)
58 write (icsa_rbank,850) (rad2deg*rvar(i,l,j,k),i=1,4)
65 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
67 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
71 end subroutine write_rbank
72 !-----------------------------------------------------------------------------
73 subroutine read_rbank(jlee,adif)
76 use geometry_data, only: nres,deg2rad
78 ! implicit real*8 (a-h,o-z)
79 ! include 'DIMENSIONS'
81 ! include 'COMMON.IOUNITS'
82 ! include 'COMMON.CSA'
83 ! include 'COMMON.BANK'
84 ! include 'COMMON.CHAIN'
85 ! include 'COMMON.GEO'
86 ! include 'COMMON.SETUP'
87 character(len=80) :: karta
89 integer :: nbankr,nstepr,nftr,icycler,kk,k,j,l,i,&
90 ierror,ierrcode,jlee,jleer
93 open(icsa_rbank,file=csa_rbank,status="old")
94 read (icsa_rbank,901) jleer,nbankr,nstepr,nftr,icycler,adif
95 print *,jleer,nbankr,nstepr,nftr,icycler,adif
96 ! print *, 'adif from read_rbank ',adif
98 if(nbankr.ne.nbank) then
99 write (iout,*) 'ERROR in READ_BANK: NBANKR',nbankr,&
101 call mpi_abort(mpi_comm_world,ierror,ierrcode)
103 if(jleer.ne.jlee) then
104 write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,&
106 call mpi_abort(mpi_comm_world,ierror,ierrcode)
112 read (icsa_rbank,'(a80)') karta
113 write(iout,*) "READ_RBANK: kk=",kk
115 ! if (index(karta,"*").gt.0) then
116 ! write (iout,*) "***** Stars in bankr ***** k=",k,
120 ! read (30,850) (rdummy,i=1,4)
125 call reada(karta,"total E",rene(kk),1.0d20)
126 call reada(karta,"rmsd from N",rrmsn(kk),0.0d0)
127 call reada(karta,"%NC",rpncn(kk),0.0d0)
128 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
129 "%NC",bpncn(kk),ibank(kk)
130 ! read (icsa_rbank,953) kdummy,rene(kk),rrmsn(kk),rpncn(kk)
133 read (icsa_rbank,850) (rvar(i,l,j,kk),i=1,4)
134 ! write (iout,850) (rvar(i,l,j,kk),i=1,4)
136 rvar(i,l,j,kk)=deg2rad*rvar(i,l,j,kk)
142 !d write (*,*) "read_rbank ******************* kk",kk,
144 if (kk.lt.nbankr) nbankr=kk
149 !d write (*,850) (rvar(i,l,j,kk),i=1,4)
156 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
157 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2)
160 end subroutine read_rbank
161 !-----------------------------------------------------------------------------
162 subroutine write_bank(jlee,nft)
165 use control_data, only: vdisulf
166 use geometry_data, only: nres,rad2deg
167 ! implicit real*8 (a-h,o-z)
168 ! include 'DIMENSIONS'
169 ! include 'COMMON.IOUNITS'
170 ! include 'COMMON.CSA'
171 ! include 'COMMON.BANK'
172 ! include 'COMMON.CHAIN'
173 ! include 'COMMON.GEO'
174 ! include 'COMMON.SBRIDGE'
175 ! include 'COMMON.CONTROL'
176 character(len=7) :: chtmp
177 character(len=40) :: chfrm
180 integer :: nft,k,l,i,j,jlee
182 open(icsa_bank,file=csa_bank,status="unknown")
183 write (icsa_bank,900) jlee,nbank,nstep,nft,icycle,cutdif
184 write (icsa_bank,902) nglob_csa, eglob_csa
185 open (igeom,file=intname,status='UNKNOWN')
187 write (icsa_bank,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
188 if (vdisulf) write (icsa_bank,'(101i4)') &
189 bvar_nss(k),((bvar_ss(j,i,k),j=1,2),i=1,bvar_nss(k))
192 write (icsa_bank,850) (rad2deg*bvar(i,l,j,k),i=1,4)
195 if (bvar_nss(k).le.9) then
196 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
197 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
199 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
200 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
201 write (igeom,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
202 bvar_ss(2,i,k),i=10,bvar_nss(k))
204 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
205 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
206 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
207 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
212 if (nstep/200.gt.ilastnstep) then
214 ilastnstep=(ilastnstep+1)*1.5
215 write(chfrm,'(a2,i1,a1)') '(i',int(dlog10(dble(nstep))+1),')'
216 write(chtmp,chfrm) nstep
217 open(icsa_int,file=prefix(:ilen(prefix)) &
218 //'_'//chtmp(:ilen(chtmp))//'.int',status='UNKNOWN')
220 if (bvar_nss(k).le.9) then
221 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
222 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,bvar_nss(k))
224 write (icsa_int,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
225 bvar_nss(k),(bvar_ss(1,i,k),bvar_ss(2,i,k),i=1,9)
226 write (icsa_int,'(3X,11(1X,2I3))') (bvar_ss(1,i,k),&
227 bvar_ss(2,i,k),i=10,bvar_nss(k))
229 write (icsa_int,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
230 write (icsa_int,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
231 write (icsa_int,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
232 write (icsa_int,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
240 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
242 902 format (1x,'nglob_csa =',i4,' eglob_csa =',1pe14.5)
243 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
247 end subroutine write_bank
248 !-----------------------------------------------------------------------------
249 subroutine write_bank_reminimized(jlee,nft)
252 use geometry_data, only: nres,rad2deg
254 ! implicit real*8 (a-h,o-z)
255 ! include 'DIMENSIONS'
256 ! include 'COMMON.IOUNITS'
257 ! include 'COMMON.CSA'
258 ! include 'COMMON.BANK'
259 ! include 'COMMON.CHAIN'
260 ! include 'COMMON.GEO'
261 ! include 'COMMON.SBRIDGE'
263 integer :: nft,i,l,j,k,jlee
265 open(icsa_bank_reminimized,file=csa_bank_reminimized,&
267 write (icsa_bank_reminimized,900) &
268 jlee,nbank,nstep,nft,icycle,cutdif
269 open (igeom,file=intname,status='UNKNOWN')
271 write (icsa_bank_reminimized,952) k,bene(k),brmsn(k),&
275 write (icsa_bank_reminimized,850) (rad2deg*bvar(i,l,j,k),i=1,4)
279 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
280 nss,(ihpb(i),jhpb(i),i=1,nss)
282 write (igeom,'(I5,F10.3,I2,9(1X,2I3))') k,bene(k),&
283 nss,(ihpb(i),jhpb(i),i=1,9)
284 write (igeom,'(3X,11(1X,2I3))') (ihpb(i),jhpb(i),i=10,nss)
286 write (igeom,200) (rad2deg*bvar(1,i,1,k),i=2,nres-1)
287 write (igeom,200) (rad2deg*bvar(2,i,1,k),i=2,nres-2)
288 write (igeom,200) (rad2deg*bvar(3,i,1,k),i=2,nres-1)
289 write (igeom,200) (rad2deg*bvar(4,i,1,k),i=2,nres-1)
291 close(icsa_bank_reminimized)
296 900 format (1x,"jlee =",i3,3x,"nbank =",i4,3x,"nstep =",&
298 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
302 end subroutine write_bank_reminimized
303 !-----------------------------------------------------------------------------
304 subroutine read_bank(jlee,nft,cutdifr)
307 use control_data, only: vdisulf
308 use geometry_data, only: nres,deg2rad
310 ! implicit real*8 (a-h,o-z)
311 ! include 'DIMENSIONS'
312 ! include 'COMMON.IOUNITS'
313 ! include 'COMMON.CSA'
314 ! include 'COMMON.BANK'
315 ! include 'COMMON.CHAIN'
316 ! include 'COMMON.GEO'
317 ! include 'COMMON.CONTROL'
318 ! include 'COMMON.SBRIDGE'
319 character(len=80) :: karta
323 integer :: nft,kk,k,l,i,j,jlee
324 real(kind=8) :: cutdifr
326 open(icsa_bank,file=csa_bank,status="old")
327 read (icsa_bank,901) jlee,nbank,nstep,nft,icycle,cutdifr
328 read (icsa_bank,902) nglob_csa, eglob_csa
329 ! if(jleer.ne.jlee) then
330 ! write (iout,*) 'ERROR in READ_BANK: JLEER',jleer,
332 ! call mpi_abort(mpi_comm_world,ierror,ierrcode)
337 read (icsa_bank,'(a80)') karta
338 write(iout,*) "READ_BANK: kk=",kk
340 ! if (index(karta,"*").gt.0) then
341 ! write (iout,*) "***** Stars in bank ***** k=",k,
345 ! read (33,850) (rdummy,i=1,4)
350 call reada(karta,"total E",bene(kk),1.0d20)
351 call reada(karta,"rmsd from N",brmsn(kk),0.0d0)
352 call reada(karta,"%NC",bpncn(kk),0.0d0)
353 read (karta(ilen(karta)-1:),*,end=111,err=111) ibank(kk)
357 write(iout,*)"total E",bene(kk),"rmsd from N",brmsn(kk),&
358 "%NC",bpncn(kk),ibank(kk)
359 ! read (icsa_bank,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
361 read (icsa_bank,'(101i4)') &
362 bvar_nss(kk),((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
363 bvar_ns(kk)=ns-2*bvar_nss(kk)
364 write(iout,*) 'read SSBOND',bvar_nss(kk),&
365 ((bvar_ss(j,i,kk),j=1,2),i=1,bvar_nss(kk))
366 !d write(iout,*) 'read CYS #free ', bvar_ns(kk)
370 do while( iss(i).ne.bvar_ss(1,j,kk)-nres .and. &
371 iss(i).ne.bvar_ss(2,j,kk)-nres .and. &
375 if (j.gt.bvar_nss(kk)) then
380 !d write(iout,*)'read CYS free',(bvar_s(l,kk),l=1,bvar_ns(kk))
384 read (icsa_bank,850) (bvar(i,l,j,kk),i=1,4)
385 ! write (iout,850) (bvar(i,l,j,kk),i=1,4)
387 bvar(i,l,j,kk)=deg2rad*bvar(i,l,j,kk)
394 if (kk.lt.nbank) nbank=kk
395 !d write (*,*) "read_bank ******************* kk",kk,
401 !d write (*,850) (bvar(i,l,j,kk),i=1,4)
407 ! read (33,953) kdummy,bene(k),brmsn(k),bpncn(k),ibank(k)
410 ! read (33,850) (bvar(i,l,j,k),i=1,4)
412 ! bvar(i,l,j,k)=deg2rad*bvar(i,l,j,k)
420 952 format (1x,'#',i4,' total E ',f12.3,' rmsd from N ',f8.3,i5)
421 901 format (1x,6x,i3,3x,7x,i4,3x,7x,i8,i10,i2,f15.5)
422 902 format (1x,11x,i4,12x,1pe14.5)
423 953 format (1x,1x,i4,9x,f12.3,13x,f8.3,5x,f5.2,i5)
426 end subroutine read_bank
427 !-----------------------------------------------------------------------------
428 subroutine write_bank1(jlee)
431 use geometry_data, only: nres,rad2deg
432 ! implicit real*8 (a-h,o-z)
433 ! include 'DIMENSIONS'
434 ! include 'COMMON.IOUNITS'
435 ! include 'COMMON.CSA'
436 ! include 'COMMON.BANK'
437 ! include 'COMMON.CHAIN'
438 ! include 'COMMON.GEO'
440 integer :: k,i,l,j,jlee
442 #if defined(AIX) || defined(PGI)
443 open(icsa_bank1,file=csa_bank1,position="append")
445 open(icsa_bank1,file=csa_bank1,access="append")
447 write (icsa_bank1,900) jlee,nbank,nstep,cutdif
449 write (icsa_bank1,952) k,bene(k),brmsn(k),bpncn(k),ibank(k)
452 write (icsa_bank1,850) (rad2deg*bvar(i,l,j,k),i=1,4)
458 900 format (4x,"jlee =",i5,3x,"nbank =",i5,3x,"nstep =",i10,f15.5)
459 952 format (1x,'#',i4,' total E ',1pe14.5,' rmsd from N ',0pf8.3,&
463 end subroutine write_bank1
464 !-----------------------------------------------------------------------------
466 !-----------------------------------------------------------------------------
467 ! subroutine cartprint
469 ! use geometry_data, only: c
470 ! use energy_data, only: itype
471 ! implicit real*8 (a-h,o-z)
472 ! include 'DIMENSIONS'
473 ! include 'COMMON.CHAIN'
474 ! include 'COMMON.INTERACT'
475 ! include 'COMMON.NAMES'
476 ! include 'COMMON.IOUNITS'
481 ! write (iout,110) restyp(itype(i,1)),i,c(1,i),c(2,i),&
482 ! c(3,i),c(1,nres+i),c(2,nres+i),c(3,nres+i)
484 ! 100 format (//' alpha-carbon coordinates ',&
485 ! ' centroid coordinates'/ &
486 ! ' ', 6X,'X',11X,'Y',11X,'Z',&
487 ! 10X,'X',11X,'Y',11X,'Z')
488 ! 110 format (a,'(',i3,')',6f12.5)
490 ! end subroutine cartprint
491 !-----------------------------------------------------------------------------
493 !-----------------------------------------------------------------------------
494 subroutine secstrp2dihc
498 ! implicit real*8 (a-h,o-z)
499 ! include 'DIMENSIONS'
500 ! include 'COMMON.GEO'
501 ! include 'COMMON.BOUNDS'
502 ! include 'COMMON.CHAIN'
503 ! include 'COMMON.TORCNSTR'
504 ! include 'COMMON.IOUNITS'
505 !el character(len=1),dimension(nres) :: secstruc !(maxres)
506 !el COMMON/SECONDARYS/secstruc
507 character(len=80) :: line
512 integer :: i,ii,lenpre
514 allocate(secstruc(nres))
516 !dr call getenv_loc('SECPREDFIL',secpred)
518 secpred=prefix(:lenpre)//'.spred'
520 #if defined(WINIFL) || defined(WINPGI)
521 open(isecpred,file=secpred,status='old',readonly,shared)
522 #elif (defined CRAY) || (defined AIX)
523 open(isecpred,file=secpred,status='old',action='read')
525 open(isecpred,file=secpred,status='old')
527 open(isecpred,file=secpred,status='old',action='read')
529 ! read secondary structure prediction from JPRED here!
530 ! read(isecpred,'(A80)',err=100,end=100) line
531 ! read(line,'(f10.3)',err=110) ftors
532 read(isecpred,'(f10.3)',err=110) ftors(1)
534 write (iout,*) 'FTORS factor =',ftors(1)
535 ! initialize secstruc to any
542 call read_secstr_pred(isecpred,iout,errflag)
544 write(iout,*)'There is a problem with the list of secondary-',&
545 'structure prediction'
548 ! 8/13/98 Set limits to generating the dihedral angles
557 if ( secstruc(i) .eq. 'H') then
558 ! Helix restraints for this residue
561 phi0(ii) = 45.0D0*deg2rad
562 drange(ii)= 5.0D0*deg2rad
563 phibound(1,i) = phi0(ii)-drange(ii)
564 phibound(2,i) = phi0(ii)+drange(ii)
565 else if (secstruc(i) .eq. 'E') then
566 ! strand restraints for this residue
569 phi0(ii) = 180.0D0*deg2rad
570 drange(ii)= 5.0D0*deg2rad
571 phibound(1,i) = phi0(ii)-drange(ii)
572 phibound(2,i) = phi0(ii)+drange(ii)
574 ! no restraints for this residue
575 ndih_nconstr=ndih_nconstr+1
576 idih_nconstr(ndih_nconstr)=i
580 ! deallocate(secstruc)
583 write(iout,'(A30,A80)')'Error reading file SECPRED',secpred
584 ! deallocate(secstruc)
587 write(iout,'(A20)')'Error reading FTORS'
588 ! deallocate(secstruc)
590 end subroutine secstrp2dihc
591 !-----------------------------------------------------------------------------
592 subroutine read_secstr_pred(jin,jout,errors)
594 ! implicit real*8 (a-h,o-z)
595 ! INCLUDE 'DIMENSIONS'
596 ! include 'COMMON.IOUNITS'
597 ! include 'COMMON.CHAIN'
598 !el character(len=1),dimension(nres) :: secstruc !(maxres)
599 !el COMMON/SECONDARYS/secstruc
601 character(len=80) :: line,line1 !,ucase
602 logical :: errflag,errors,blankline
605 integer :: jin,jout,iseq,ipos,ipos1,iend,il,&
608 read (jin,'(a)') line
609 write (jout,'(2a)') '> ',line(1:78)
611 ! Remember that we number full residues starting from 2, then, iseq=1 and iseq=nres
612 ! correspond to the end-groups. ADD to the secondary structure prediction "-" for the
613 ! end-groups in the input file "*.spred"
616 do while (index(line1,'$END').eq.0)
617 ! Override commented lines.
620 do while (.not.blankline)
622 call mykey(line,line1,ipos,blankline,errflag)
623 if (errflag) write (jout,'(2a)') &
624 'Error when reading sequence in line: ',line
625 errors=errors .or. errflag
626 if (.not. blankline .and. .not. errflag) then
629 !el if (iseq.le.maxres) then
630 if (line1(1:1).eq.'-' ) then
631 secstruc(iseq)=line1(1:1)
632 else if ( ( ucase(line1(1:1)).eq.'E' ) .or. &
633 ( ucase(line1(1:1)).eq.'H' ) ) then
634 secstruc(iseq)=ucase(line1(1:1))
637 write (jout,1010) line1(1:1), iseq
642 !el write (jout,1000) iseq,maxres
645 do while (ipos1.le.iend)
650 !el if (iseq.le.maxres) then
651 if (line1(ipos1-1:ipos1-1).eq.'-' ) then
652 secstruc(iseq)=line1(ipos1-1:ipos1-1)
653 else if((ucase(line1(ipos1-1:ipos1-1)).eq.'E').or. &
654 (ucase(line1(ipos1-1:ipos1-1)).eq.'H') ) then
655 secstruc(iseq)=ucase(line1(ipos1-1:ipos1-1))
658 write (jout,1010) line1(ipos1-1:ipos1-1), iseq
663 !el write (jout,1000) iseq,maxres
670 read (jin,'(a)') line
671 write (jout,'(2a)') '> ',line(1:78)
675 !d write (jout,'(10a8)') (sequence(i),i=1,iseq-1)
677 !d check whether the found length of the chain is correct.
678 length_of_chain=iseq-1
679 if (length_of_chain .ne. nres) then
681 write (jout,'(a,i4,a,i4,a)') &
682 'Error: the number of labels specified in $SEC_STRUC_PRED (' &
683 ,length_of_chain,') does not match with the number of residues (' &
688 1000 format('Error - the number of residues (',i4,&
689 ') has exceeded maximum (',i4,').')
690 1010 format ('Error - unrecognized secondary structure label',a4,&
693 end subroutine read_secstr_pred
695 !-----------------------------------------------------------------------------
697 !-----------------------------------------------------------------------------
702 use control_data, only:maxterm !,maxtor
706 use control, only: getenv_loc
708 ! Read the parameters of the probability distributions of the virtual-bond
709 ! valence angles and the side chains and energy parameters.
711 ! Important! Energy-term weights ARE NOT read here; they are read from the
712 ! main input file instead, because NO defaults have yet been set for these
715 ! implicit real*8 (a-h,o-z)
716 ! include 'DIMENSIONS'
721 ! include 'COMMON.IOUNITS'
722 ! include 'COMMON.CHAIN'
723 ! include 'COMMON.INTERACT'
724 ! include 'COMMON.GEO'
725 ! include 'COMMON.LOCAL'
726 ! include 'COMMON.TORSION'
727 ! include 'COMMON.SCCOR'
728 ! include 'COMMON.SCROT'
729 ! include 'COMMON.FFIELD'
730 ! include 'COMMON.NAMES'
731 ! include 'COMMON.SBRIDGE'
732 ! include 'COMMON.MD'
733 ! include 'COMMON.SETUP'
734 character(len=1) :: t1,t2,t3
735 character(len=1) :: onelett(4) = (/"G","A","P","D"/)
736 character(len=1) :: toronelet(-2:2) = (/"p","a","G","A","P"/)
737 logical :: lprint,LaTeX,SPLIT_FOURIERTOR
738 real(kind=8),dimension(3,3,maxlob) :: blower !(3,3,maxlob)
739 real(kind=8),dimension(13) :: buse
740 character(len=3) :: lancuch !,ucase
742 integer :: m,n,l,i,j,k,iblock,lll,llll,ll,nlobi,mm,jj
743 integer :: maxinter,junk,kk,ii,ncatprotparm,nkcctyp
744 real(kind=8) :: v0ijsccor,v0ijsccor1,v0ijsccor2,v0ijsccor3,si,&
745 dwa16,rjunk,akl,v0ij,rri,epsij,rrij,sigeps,sigt1sq,&
746 sigt2sq,sigii1,sigii2,ratsig1,ratsig2,rsum_max,r_augm,&
747 res1,epsijlip,epspeptube,epssctube,sigmapeptube, &
749 integer :: ichir1,ichir2,ijunk
752 ! real(kind=8),dimension(maxterm,-maxtor:maxtor,-maxtor:maxtor,2) :: v1_el,v2_el !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
753 !el allocate(v1_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
754 !el allocate(v2_el(maxterm,-maxtor:maxtor,-maxtor:maxtor,2))
756 ! For printing parameters after they are read set the following in the UNRES
759 ! setenv PRINT_PARM YES
761 ! To print parameters in LaTeX format rather than as ASCII tables:
765 call getenv_loc("PRINT_PARM",lancuch)
766 lprint = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
767 call getenv_loc("LATEX",lancuch)
768 LaTeX = (ucase(lancuch).eq."YES" .or. ucase(lancuch).eq."Y")
770 dwa16=2.0d0**(1.0d0/6.0d0)
772 ! Assign virtual-bond length
775 vblinv2=vblinv*vblinv
777 ! Read the virtual-bond parameters, masses, and moments of inertia
778 ! and Stokes' radii of the peptide group and side chains
780 allocate(dsc(ntyp1)) !(ntyp1)
781 allocate(dsc_inv(ntyp1)) !(ntyp1)
782 allocate(nbondterm_nucl(ntyp_molec(2))) !(ntyp)
783 allocate(vbldsc0_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
784 allocate(aksc_nucl(maxbondterm,ntyp_molec(2))) !(maxbondterm,ntyp)
785 allocate(nbondterm(ntyp)) !(ntyp)
786 allocate(vbldsc0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
787 allocate(aksc(maxbondterm,ntyp)) !(maxbondterm,ntyp)
788 allocate(abond0(maxbondterm,ntyp)) !(maxbondterm,ntyp)
789 allocate(long_r_sidechain(ntyp))
790 allocate(short_r_sidechain(ntyp))
795 allocate(msc(ntyp+1)) !(ntyp+1)
796 allocate(isc(ntyp+1)) !(ntyp+1)
797 allocate(restok(ntyp+1)) !(ntyp+1)
799 read (ibond,*) vbldp0,akp,mp,ip,pstok
802 read (ibond,*) vbldsc0(1,i),aksc(1,i),msc(i),isc(i),restok(i)
803 dsc(i) = vbldsc0(1,i)
807 dsc_inv(i)=1.0D0/dsc(i)
816 allocate(msc(ntyp+1,5)) !(ntyp+1)
817 allocate(isc(ntyp+1,5)) !(ntyp+1)
818 allocate(restok(ntyp+1,5)) !(ntyp+1)
820 read (ibond,*) junk,vbldp0,vbldpDUM,akp,rjunk,mp(1),ip(1),pstok(1)
822 read (ibond,*) nbondterm(i),(vbldsc0(j,i),aksc(j,i),abond0(j,i),&
823 j=1,nbondterm(i)),msc(i,1),isc(i,1),restok(i,1)
824 dsc(i) = vbldsc0(1,i)
828 dsc_inv(i)=1.0D0/dsc(i)
832 read (ibond_nucl,*) vbldp0_nucl,akp_nucl,mp(2),ip(2),pstok(2)
835 read (ibond_nucl,*) vbldsc0_nucl(1,i),aksc_nucl(1,i),msc(i,2),isc(i,2),restok(i,2)
836 ! dsc(i) = vbldsc0_nucl(1,i)
840 ! dsc_inv(i)=1.0D0/dsc(i)
843 ! read (ibond_nucl,*) junk,vbldp0_nucl,akp_nucl,rjunk,mp(2),ip(2),pstok(2)
844 ! do i=1,ntyp_molec(2)
845 ! read (ibond_nucl,*) nbondterm_nucl(i),(vbldsc0_nucl(j,i),&
846 ! aksc_nucl(j,i),abond0_nucl(j,i),&
847 ! j=1,nbondterm_nucl(i)),msc(i,2),isc(i,2),restok(i,2)
848 ! dsc(i) = vbldsc0(1,i)
852 ! dsc_inv(i)=1.0D0/dsc(i)
857 write(iout,'(/a/)')"Dynamic constants of the interaction sites:"
858 write (iout,'(a10,a3,6a10)') 'Type','N','VBL','K','A0','mass',&
860 write(iout,'(a10,i3,6f10.5)') "p",1,vbldp0,akp,0.0d0,mp(1),ip(1),pstok(1)
862 write (iout,'(a10,i3,6f10.5)') restyp(i,1),nbondterm(i),&
863 vbldsc0(1,i),aksc(1,i),abond0(1,i),msc(i,1),isc(i,1),restok(i,1)
865 write (iout,'(13x,3f10.5)') &
866 vbldsc0(j,i),aksc(j,i),abond0(j,i)
871 read(iion,*) msc(i,5),restok(i,5)
872 print *,msc(i,5),restok(i,5)
876 read (iion,*) ncatprotparm
877 allocate(catprm(ncatprotparm,4))
879 read (iion,*) (catprm(i,k),i=1,ncatprotparm)
882 ! read (iion,*) (vcatprm(k),k=1,ncatprotpram)
883 !----------------------------------------------------
884 allocate(a0thet(-ntyp:ntyp),theta0(-ntyp:ntyp))
885 allocate(sig0(-ntyp:ntyp),sigc0(-ntyp:ntyp)) !(-ntyp:ntyp)
886 allocate(athet(2,-ntyp:ntyp,-1:1,-1:1))
887 allocate(bthet(2,-ntyp:ntyp,-1:1,-1:1)) !(2,-ntyp:ntyp,-1:1,-1:1)
888 allocate(polthet(0:3,-ntyp:ntyp)) !(0:3,-ntyp:ntyp)
889 allocate(gthet(3,-ntyp:ntyp)) !(3,-ntyp:ntyp)
899 allocate(liptranene(ntyp))
900 !C reading lipid parameters
901 write (iout,*) "iliptranpar",iliptranpar
903 read(iliptranpar,*) pepliptran
906 read(iliptranpar,*) liptranene(i)
907 print *,liptranene(i)
913 ! Read the parameters of the probability distribution/energy expression
914 ! of the virtual-bond valence angles theta
917 read (ithep,*,err=111,end=111) a0thet(i),(athet(j,i,1,1),j=1,2),&
918 (bthet(j,i,1,1),j=1,2)
919 read (ithep,*,err=111,end=111) (polthet(j,i),j=0,3)
920 read (ithep,*,err=111,end=111) (gthet(j,i),j=1,3)
921 read (ithep,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
925 athet(1,i,1,-1)=athet(1,i,1,1)
926 athet(2,i,1,-1)=athet(2,i,1,1)
927 bthet(1,i,1,-1)=-bthet(1,i,1,1)
928 bthet(2,i,1,-1)=-bthet(2,i,1,1)
929 athet(1,i,-1,1)=-athet(1,i,1,1)
930 athet(2,i,-1,1)=-athet(2,i,1,1)
931 bthet(1,i,-1,1)=bthet(1,i,1,1)
932 bthet(2,i,-1,1)=bthet(2,i,1,1)
936 athet(1,i,-1,-1)=athet(1,-i,1,1)
937 athet(2,i,-1,-1)=-athet(2,-i,1,1)
938 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
939 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
940 athet(1,i,-1,1)=athet(1,-i,1,1)
941 athet(2,i,-1,1)=-athet(2,-i,1,1)
942 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
943 bthet(2,i,-1,1)=bthet(2,-i,1,1)
944 athet(1,i,1,-1)=-athet(1,-i,1,1)
945 athet(2,i,1,-1)=athet(2,-i,1,1)
946 bthet(1,i,1,-1)=bthet(1,-i,1,1)
947 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
952 polthet(j,i)=polthet(j,-i)
955 gthet(j,i)=gthet(j,-i)
963 'Parameters of the virtual-bond valence angles:'
964 write (iout,'(/a/9x,5a/79(1h-))') 'Fourier coefficients:',&
965 ' ATHETA0 ',' A1 ',' A2 ',&
968 write(iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
969 a0thet(i),(athet(j,i,1,1),j=1,2),(bthet(j,i,1,1),j=1,2)
971 write (iout,'(/a/9x,5a/79(1h-))') &
972 'Parameters of the expression for sigma(theta_c):',&
973 ' ALPH0 ',' ALPH1 ',' ALPH2 ',&
974 ' ALPH3 ',' SIGMA0C '
976 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,&
977 (polthet(j,i),j=0,3),sigc0(i)
979 write (iout,'(/a/9x,5a/79(1h-))') &
980 'Parameters of the second gaussian:',&
981 ' THETA0 ',' SIGMA0 ',' G1 ',&
984 write (iout,'(a3,i4,2x,5(1pe14.5))') restyp(i,1),i,theta0(i),&
985 sig0(i),(gthet(j,i),j=1,3)
989 'Parameters of the virtual-bond valence angles:'
990 write (iout,'(/a/9x,5a/79(1h-))') &
991 'Coefficients of expansion',&
992 ' theta0 ',' a1*10^2 ',' a2*10^2 ',&
993 ' b1*10^1 ',' b2*10^1 '
995 write(iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),&
996 a0thet(i),(100*athet(j,i,1,1),j=1,2),&
997 (10*bthet(j,i,1,1),j=1,2)
999 write (iout,'(/a/9x,5a/79(1h-))') &
1000 'Parameters of the expression for sigma(theta_c):',&
1001 ' alpha0 ',' alph1 ',' alph2 ',&
1002 ' alhp3 ',' sigma0c '
1004 write (iout,'(a3,1h&,2x,5(1pe12.3,1h&))') restyp(i,1),&
1005 (polthet(j,i),j=0,3),sigc0(i)
1007 write (iout,'(/a/9x,5a/79(1h-))') &
1008 'Parameters of the second gaussian:',&
1009 ' theta0 ',' sigma0*10^2 ',' G1*10^-1',&
1012 write (iout,'(a3,1h&,2x,5(f8.3,1h&))') restyp(i,1),theta0(i),&
1013 100*sig0(i),gthet(1,i)*0.1D0,gthet(2,i),gthet(3,i)*10.0D0
1019 ! Read the parameters of Utheta determined from ab initio surfaces
1020 ! Kozlowska et al., J. Phys.: Condens. Matter 19 (2007) 285203
1022 IF (tor_mode.eq.0) THEN
1023 read (ithep,*,err=111,end=111) nthetyp,ntheterm,ntheterm2,&
1024 ntheterm3,nsingle,ndouble
1025 nntheterm=max0(ntheterm,ntheterm2,ntheterm3)
1027 !----------------------------------------------------
1028 allocate(ithetyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
1029 allocate(aa0thet(-nthetyp-1:nthetyp+1,&
1030 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1031 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1032 allocate(aathet(ntheterm,-nthetyp-1:nthetyp+1,&
1033 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1034 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1035 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1036 allocate(bbthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1037 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1038 allocate(ccthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1039 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1040 allocate(ddthet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1041 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1042 allocate(eethet(nsingle,ntheterm2,-nthetyp-1:nthetyp+1,&
1043 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1044 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1045 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1046 allocate(ffthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1047 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1048 allocate(ggthet(ndouble,ndouble,ntheterm3,-nthetyp-1:nthetyp+1,&
1049 -nthetyp-1:nthetyp+1,-nthetyp-1:nthetyp+1,2))
1050 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1051 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1053 read (ithep,*,err=111,end=111) (ithetyp(i),i=1,ntyp1)
1055 ithetyp(i)=-ithetyp(-i)
1058 aa0thet(:,:,:,:)=0.0d0
1059 aathet(:,:,:,:,:)=0.0d0
1060 bbthet(:,:,:,:,:,:)=0.0d0
1061 ccthet(:,:,:,:,:,:)=0.0d0
1062 ddthet(:,:,:,:,:,:)=0.0d0
1063 eethet(:,:,:,:,:,:)=0.0d0
1064 ffthet(:,:,:,:,:,:,:)=0.0d0
1065 ggthet(:,:,:,:,:,:,:)=0.0d0
1067 ! VAR:iblock means terminally blocking group 1=non-proline 2=proline
1069 ! VAR:ntethtyp is type of theta potentials type currently 0=glycine
1070 ! VAR:1=non-glicyne non-proline 2=proline
1071 ! VAR:negative values for D-aminoacid
1073 do j=-nthetyp,nthetyp
1074 do k=-nthetyp,nthetyp
1075 read (ithep,'(6a)',end=111,err=111) res1
1076 read (ithep,*,end=111,err=111) aa0thet(i,j,k,iblock)
1077 ! VAR: aa0thet is variable describing the average value of Foureir
1078 ! VAR: expansion series
1079 ! VAR: aathet is foureir expansion in theta/2 angle for full formula
1080 ! VAR: look at the fitting equation in Kozlowska et al., J. Phys.:
1081 !ondens. Matter 19 (2007) 285203 and Sieradzan et al., unpublished
1082 read (ithep,*,end=111,err=111) &
1083 (aathet(l,i,j,k,iblock),l=1,ntheterm)
1084 read (ithep,*,end=111,err=111) &
1085 ((bbthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1086 (ccthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1087 (ddthet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1088 (eethet(lll,ll,i,j,k,iblock),lll=1,nsingle),&
1090 read (ithep,*,end=111,err=111) &
1091 (((ffthet(llll,lll,ll,i,j,k,iblock),&
1092 ffthet(lll,llll,ll,i,j,k,iblock),&
1093 ggthet(llll,lll,ll,i,j,k,iblock),&
1094 ggthet(lll,llll,ll,i,j,k,iblock),&
1095 llll=1,lll-1),lll=2,ndouble),ll=1,ntheterm3)
1100 ! For dummy ends assign glycine-type coefficients of theta-only terms; the
1101 ! coefficients of theta-and-gamma-dependent terms are zero.
1102 ! IF YOU WANT VALENCE POTENTIALS FOR DUMMY ATOM UNCOMENT BELOW (NOT
1103 ! RECOMENTDED AFTER VERSION 3.3)
1107 ! aathet(l,i,j,nthetyp+1,iblock)=aathet(l,i,j,1,iblock)
1108 ! aathet(l,nthetyp+1,i,j,iblock)=aathet(l,1,i,j,iblock)
1110 ! aa0thet(i,j,nthetyp+1,iblock)=aa0thet(i,j,1,iblock)
1111 ! aa0thet(nthetyp+1,i,j,iblock)=aa0thet(1,i,j,iblock)
1114 ! aathet(l,nthetyp+1,i,nthetyp+1,iblock)=aathet(l,1,i,1,iblock)
1116 ! aa0thet(nthetyp+1,i,nthetyp+1,iblock)=aa0thet(1,i,1,iblock)
1119 ! AND COMMENT THE LOOPS BELOW
1123 aathet(l,i,j,nthetyp+1,iblock)=0.0d0
1124 aathet(l,nthetyp+1,i,j,iblock)=0.0d0
1126 aa0thet(i,j,nthetyp+1,iblock)=0.0d0
1127 aa0thet(nthetyp+1,i,j,iblock)=0.0d0
1130 aathet(l,nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1132 aa0thet(nthetyp+1,i,nthetyp+1,iblock)=0.0d0
1137 ! Substitution for D aminoacids from symmetry.
1140 do j=-nthetyp,nthetyp
1141 do k=-nthetyp,nthetyp
1142 aa0thet(i,j,k,iblock)=aa0thet(-i,-j,-k,iblock)
1144 aathet(l,i,j,k,iblock)=aathet(l,-i,-j,-k,iblock)
1148 bbthet(lll,ll,i,j,k,iblock)=bbthet(lll,ll,-i,-j,-k,iblock)
1149 ccthet(lll,ll,i,j,k,iblock)=-ccthet(lll,ll,-i,-j,-k,iblock)
1150 ddthet(lll,ll,i,j,k,iblock)=ddthet(lll,ll,-i,-j,-k,iblock)
1151 eethet(lll,ll,i,j,k,iblock)=-eethet(lll,ll,-i,-j,-k,iblock)
1157 ffthet(llll,lll,ll,i,j,k,iblock)= &
1158 ffthet(llll,lll,ll,-i,-j,-k,iblock)
1159 ffthet(lll,llll,ll,i,j,k,iblock)= &
1160 ffthet(lll,llll,ll,-i,-j,-k,iblock)
1161 ggthet(llll,lll,ll,i,j,k,iblock)= &
1162 -ggthet(llll,lll,ll,-i,-j,-k,iblock)
1163 ggthet(lll,llll,ll,i,j,k,iblock)= &
1164 -ggthet(lll,llll,ll,-i,-j,-k,iblock)
1173 ! Control printout of the coefficients of virtual-bond-angle potentials
1176 write (iout,'(//a)') 'Parameter of virtual-bond-angle potential'
1181 write (iout,'(//4a)') &
1182 'Type ',onelett(i),onelett(j),onelett(k)
1183 write (iout,'(//a,10x,a)') " l","a[l]"
1184 write (iout,'(i2,1pe15.5)') 0,aa0thet(i,j,k,iblock)
1185 write (iout,'(i2,1pe15.5)') &
1186 (l,aathet(l,i,j,k,iblock),l=1,ntheterm)
1188 write (iout,'(//2h m,4(9x,a,3h[m,,i1,1h]))') &
1189 "b",l,"c",l,"d",l,"e",l
1191 write (iout,'(i2,4(1pe15.5))') m,&
1192 bbthet(m,l,i,j,k,iblock),ccthet(m,l,i,j,k,iblock),&
1193 ddthet(m,l,i,j,k,iblock),eethet(m,l,i,j,k,iblock)
1197 write (iout,'(//3hm,n,4(6x,a,5h[m,n,,i1,1h]))') &
1198 "f+",l,"f-",l,"g+",l,"g-",l
1201 write (iout,'(i1,1x,i1,4(1pe15.5))') n,m,&
1202 ffthet(n,m,l,i,j,k,iblock),&
1203 ffthet(m,n,l,i,j,k,iblock),&
1204 ggthet(n,m,l,i,j,k,iblock),&
1205 ggthet(m,n,l,i,j,k,iblock)
1216 !C here will be the apropriate recalibrating for D-aminoacid
1217 read (ithep,*,end=121,err=121) nthetyp
1218 allocate(nbend_kcc_Tb(-nthetyp:nthetyp))
1219 allocate(v1bend_chyb(0:36,-nthetyp:nthetyp))
1220 do i=-nthetyp+1,nthetyp-1
1221 read (ithep,*,end=121,err=121) nbend_kcc_Tb(i)
1222 do j=0,nbend_kcc_Tb(i)
1223 read (ithep,*,end=121,err=121) ijunk,v1bend_chyb(j,i)
1227 write (iout,'(a)') &
1228 "Parameters of the valence-only potentials"
1229 do i=-nthetyp+1,nthetyp-1
1230 write (iout,'(2a)') "Type ",toronelet(i)
1231 do k=0,nbend_kcc_Tb(i)
1232 write(iout,'(i5,f15.5)') k,v1bend_chyb(k,i)
1238 write (2,*) "Start reading THETA_PDB",ithep_pdb
1240 ! write (2,*) 'i=',i
1241 read (ithep_pdb,*,err=111,end=111) &
1242 a0thet(i),(athet(j,i,1,1),j=1,2),&
1243 (bthet(j,i,1,1),j=1,2)
1244 read (ithep_pdb,*,err=111,end=111) (polthet(j,i),j=0,3)
1245 read (ithep_pdb,*,err=111,end=111) (gthet(j,i),j=1,3)
1246 read (ithep_pdb,*,err=111,end=111) theta0(i),sig0(i),sigc0(i)
1247 sigc0(i)=sigc0(i)**2
1250 athet(1,i,1,-1)=athet(1,i,1,1)
1251 athet(2,i,1,-1)=athet(2,i,1,1)
1252 bthet(1,i,1,-1)=-bthet(1,i,1,1)
1253 bthet(2,i,1,-1)=-bthet(2,i,1,1)
1254 athet(1,i,-1,1)=-athet(1,i,1,1)
1255 athet(2,i,-1,1)=-athet(2,i,1,1)
1256 bthet(1,i,-1,1)=bthet(1,i,1,1)
1257 bthet(2,i,-1,1)=bthet(2,i,1,1)
1260 a0thet(i)=a0thet(-i)
1261 athet(1,i,-1,-1)=athet(1,-i,1,1)
1262 athet(2,i,-1,-1)=-athet(2,-i,1,1)
1263 bthet(1,i,-1,-1)=bthet(1,-i,1,1)
1264 bthet(2,i,-1,-1)=-bthet(2,-i,1,1)
1265 athet(1,i,-1,1)=athet(1,-i,1,1)
1266 athet(2,i,-1,1)=-athet(2,-i,1,1)
1267 bthet(1,i,-1,1)=-bthet(1,-i,1,1)
1268 bthet(2,i,-1,1)=bthet(2,-i,1,1)
1269 athet(1,i,1,-1)=-athet(1,-i,1,1)
1270 athet(2,i,1,-1)=athet(2,-i,1,1)
1271 bthet(1,i,1,-1)=bthet(1,-i,1,1)
1272 bthet(2,i,1,-1)=-bthet(2,-i,1,1)
1273 theta0(i)=theta0(-i)
1277 polthet(j,i)=polthet(j,-i)
1280 gthet(j,i)=gthet(j,-i)
1283 write (2,*) "End reading THETA_PDB"
1287 !--------------- Reading theta parameters for nucleic acid-------
1288 read (ithep_nucl,*,err=111,end=111) nthetyp_nucl,ntheterm_nucl,&
1289 ntheterm2_nucl,ntheterm3_nucl,nsingle_nucl,ndouble_nucl
1290 nntheterm_nucl=max0(ntheterm_nucl,ntheterm2_nucl,ntheterm3_nucl)
1291 allocate(ithetyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
1292 allocate(aa0thet_nucl(nthetyp_nucl+1,&
1293 nthetyp_nucl+1,nthetyp_nucl+1))
1294 !(-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1295 allocate(aathet_nucl(ntheterm_nucl+1,nthetyp_nucl+1,&
1296 nthetyp_nucl+1,nthetyp_nucl+1))
1297 !(maxtheterm,-maxthetyp1:maxthetyp1,&
1298 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1299 allocate(bbthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1300 nthetyp_nucl+1,nthetyp_nucl+1))
1301 allocate(ccthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1302 nthetyp_nucl+1,nthetyp_nucl+1))
1303 allocate(ddthet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1304 nthetyp_nucl+1,nthetyp_nucl+1))
1305 allocate(eethet_nucl(nsingle_nucl+1,ntheterm2_nucl+1,nthetyp_nucl+1,&
1306 nthetyp_nucl+1,nthetyp_nucl+1))
1307 !(maxsingle,maxtheterm2,-maxthetyp1:maxthetyp1,&
1308 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2)
1309 allocate(ffthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1310 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1311 allocate(ggthet_nucl(ndouble_nucl+1,ndouble_nucl+1,ntheterm3_nucl+1,&
1312 nthetyp_nucl+1,nthetyp_nucl+1,nthetyp_nucl+1))
1314 !(maxdouble,maxdouble,maxtheterm3,-maxthetyp1:maxthetyp1,&
1315 ! -maxthetyp1:maxthetyp1,-maxthetyp1:maxthetyp1,2))
1317 read (ithep_nucl,*,err=111,end=111) (ithetyp_nucl(i),i=1,ntyp1_molec(2))
1319 aa0thet_nucl(:,:,:)=0.0d0
1320 aathet_nucl(:,:,:,:)=0.0d0
1321 bbthet_nucl(:,:,:,:,:)=0.0d0
1322 ccthet_nucl(:,:,:,:,:)=0.0d0
1323 ddthet_nucl(:,:,:,:,:)=0.0d0
1324 eethet_nucl(:,:,:,:,:)=0.0d0
1325 ffthet_nucl(:,:,:,:,:,:)=0.0d0
1326 ggthet_nucl(:,:,:,:,:,:)=0.0d0
1331 read (ithep_nucl,'(3a)',end=111,err=111) t1,t2,t3
1332 read (ithep_nucl,*,end=111,err=111) aa0thet_nucl(i,j,k)
1333 read (ithep_nucl,*,end=111,err=111)(aathet_nucl(l,i,j,k),l=1,ntheterm_nucl)
1334 read (ithep_nucl,*,end=111,err=111) &
1335 (((bbthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1336 (ccthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1337 (ddthet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl), &
1338 (eethet_nucl(lll,ll,i,j,k),lll=1,nsingle_nucl)),ll=1,ntheterm2_nucl)
1339 read (ithep_nucl,*,end=111,err=111) &
1340 (((ffthet_nucl(llll,lll,ll,i,j,k),ffthet_nucl(lll,llll,ll,i,j,k), &
1341 ggthet_nucl(llll,lll,ll,i,j,k),ggthet_nucl(lll,llll,ll,i,j,k), &
1342 llll=1,lll-1),lll=2,ndouble_nucl),ll=1,ntheterm3_nucl)
1347 !-------------------------------------------
1348 allocate(nlob(ntyp1)) !(ntyp1)
1349 allocate(bsc(maxlob,ntyp)) !(maxlob,ntyp)
1350 allocate(censc(3,maxlob,-ntyp:ntyp)) !(3,maxlob,-ntyp:ntyp)
1351 allocate(gaussc(3,3,maxlob,-ntyp:ntyp)) !(3,3,maxlob,-ntyp:ntyp)
1358 gaussc(:,:,:,:)=0.0D0
1362 ! Read the parameters of the probability distribution/energy expression
1363 ! of the side chains.
1366 read (irotam,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1370 dsc_inv(i)=1.0D0/dsc(i)
1381 read(irotam,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1382 ((blower(k,l,1),l=1,k),k=1,3)
1383 censc(1,1,-i)=censc(1,1,i)
1384 censc(2,1,-i)=censc(2,1,i)
1385 censc(3,1,-i)=-censc(3,1,i)
1387 read (irotam,*,end=112,err=112) bsc(j,i)
1388 read (irotam,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1389 ((blower(k,l,j),l=1,k),k=1,3)
1390 censc(1,j,-i)=censc(1,j,i)
1391 censc(2,j,-i)=censc(2,j,i)
1392 censc(3,j,-i)=-censc(3,j,i)
1393 ! BSC is amplitude of Gaussian
1400 akl=akl+blower(k,m,j)*blower(l,m,j)
1404 if (((k.eq.3).and.(l.ne.3)) &
1405 .or.((l.eq.3).and.(k.ne.3))) then
1406 gaussc(k,l,j,-i)=-akl
1407 gaussc(l,k,j,-i)=-akl
1409 gaussc(k,l,j,-i)=akl
1410 gaussc(l,k,j,-i)=akl
1419 write (iout,'(/a)') 'Parameters of side-chain local geometry'
1422 if (nlobi.gt.0) then
1424 write (iout,'(/3a,i2,a,f8.3)') 'Residue type: ',restyp(i,1),&
1425 ' # of gaussian lobes:',nlobi,' dsc:',dsc(i)
1426 write (iout,'(1h&,a,3(2h&&,f8.3,2h&&))') &
1427 'log h',(bsc(j,i),j=1,nlobi)
1428 write (iout,'(1h&,a,3(1h&,f8.3,1h&,f8.3,1h&,f8.3,1h&))') &
1429 'x',((censc(k,j,i),k=1,3),j=1,nlobi)
1431 write (iout,'(2h& ,5(2x,1h&,3(f7.3,1h&)))') &
1432 ((gaussc(k,l,j,i),l=1,3),j=1,nlobi)
1435 write (iout,'(/a,8x,i1,4(25x,i1))') 'Lobe:',(j,j=1,nlobi)
1436 write (iout,'(a,f10.4,4(16x,f10.4))') &
1437 'Center ',(bsc(j,i),j=1,nlobi)
1438 write (iout,'(5(2x,3f8.4))') ((censc(k,j,i),k=1,3),&
1447 ! Read scrot parameters for potentials determined from all-atom AM1 calculations
1448 ! added by Urszula Kozlowska 07/11/2007
1450 !el Maximum number of SC local term fitting function coefficiants
1451 !el integer,parameter :: maxsccoef=65
1453 allocate(sc_parmin(65,ntyp)) !(maxsccoef,ntyp)
1456 read (irotam,*,end=112,err=112)
1458 read (irotam,*,end=112,err=112)
1461 read(irotam,*,end=112,err=112) sc_parmin(j,i)
1465 !---------reading nucleic acid parameters for rotamers-------------------
1466 allocate(sc_parmin_nucl(9,ntyp_molec(2))) !(maxsccoef,ntyp)
1467 do i=1,ntyp_molec(2)
1468 read (irotam_nucl,*,end=112,err=112)
1470 read(irotam_nucl,*,end=112,err=112) sc_parmin_nucl(j,i)
1476 write (iout,*) "Base rotamer parameters"
1477 do i=1,ntyp_molec(2)
1478 write (iout,'(a)') restyp(i,2)
1479 write (iout,'(i5,f10.5)') (i,sc_parmin_nucl(j,i),j=1,9)
1484 ! Read the parameters of the probability distribution/energy expression
1485 ! of the side chains.
1487 write (2,*) "Start reading ROTAM_PDB"
1489 read (irotam_pdb,'(3x,i3,f8.3)',end=112,err=112) nlob(i),dsc(i)
1493 dsc_inv(i)=1.0D0/dsc(i)
1504 read(irotam_pdb,*,end=112,err=112)(censc(k,1,i),k=1,3),&
1505 ((blower(k,l,1),l=1,k),k=1,3)
1507 read (irotam_pdb,*,end=112,err=112) bsc(j,i)
1508 read (irotam_pdb,*,end=112,err=112) (censc(k,j,i),k=1,3),&
1509 ((blower(k,l,j),l=1,k),k=1,3)
1516 akl=akl+blower(k,m,j)*blower(l,m,j)
1526 write (2,*) "End reading ROTAM_PDB"
1532 !C 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
1533 !C interaction energy of the Gly, Ala, and Pro prototypes.
1535 read (ifourier,*) nloctyp
1536 SPLIT_FOURIERTOR = nloctyp.lt.0
1537 nloctyp = iabs(nloctyp)
1538 !C allocate(b1(2,nres)) !(2,-maxtor:maxtor)
1539 !C allocate(b2(2,nres)) !(2,-maxtor:maxtor)
1540 !C allocate(b1tilde(2,nres)) !(2,-maxtor:maxtor)
1541 !C allocate(ctilde(2,2,nres))
1542 !C allocate(dtilde(2,2,nres)) !(2,2,-maxtor:maxtor)
1543 !C allocate(gtb1(2,nres))
1544 !C allocate(gtb2(2,nres))
1545 !C allocate(cc(2,2,nres))
1546 !C allocate(dd(2,2,nres))
1547 !C allocate(ee(2,2,nres))
1548 !C allocate(gtcc(2,2,nres))
1549 !C allocate(gtdd(2,2,nres))
1550 !C allocate(gtee(2,2,nres))
1553 allocate(itype2loc(-ntyp1:ntyp1))
1554 allocate(iloctyp(-nloctyp:nloctyp))
1555 allocate(bnew1(3,2,-nloctyp:nloctyp))
1556 allocate(bnew2(3,2,-nloctyp:nloctyp))
1557 allocate(ccnew(3,2,-nloctyp:nloctyp))
1558 allocate(ddnew(3,2,-nloctyp:nloctyp))
1559 allocate(e0new(3,-nloctyp:nloctyp))
1560 allocate(eenew(2,2,2,-nloctyp:nloctyp))
1561 allocate(bnew1tor(3,2,-nloctyp:nloctyp))
1562 allocate(bnew2tor(3,2,-nloctyp:nloctyp))
1563 allocate(ccnewtor(3,2,-nloctyp:nloctyp))
1564 allocate(ddnewtor(3,2,-nloctyp:nloctyp))
1565 allocate(e0newtor(3,-nloctyp:nloctyp))
1566 allocate(eenewtor(2,2,2,-nloctyp:nloctyp))
1568 read (ifourier,*,end=115,err=115) (itype2loc(i),i=1,ntyp)
1569 read (ifourier,*,end=115,err=115) (iloctyp(i),i=0,nloctyp-1)
1570 itype2loc(ntyp1)=nloctyp
1571 iloctyp(nloctyp)=ntyp1
1573 itype2loc(-i)=-itype2loc(i)
1576 allocate(iloctyp(-nloctyp:nloctyp))
1577 allocate(itype2loc(-ntyp1:ntyp1))
1584 iloctyp(-i)=-iloctyp(i)
1586 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1587 !c write (iout,*) "nloctyp",nloctyp,
1588 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1589 !c write (iout,*) "itype2loc",(itype2loc(i),i=1,ntyp1)
1590 !c write (iout,*) "nloctyp",nloctyp,
1591 !c & " iloctyp",(iloctyp(i),i=0,nloctyp)
1594 !c write (iout,*) "NEWCORR",i
1595 read (ifourier,*,end=115,err=115)
1598 read (ifourier,*,end=115,err=115) bnew1(ii,j,i)
1601 !c write (iout,*) "NEWCORR BNEW1"
1602 !c write (iout,*) ((bnew1(ii,j,i),ii=1,3),j=1,2)
1605 read (ifourier,*,end=115,err=115) bnew2(ii,j,i)
1608 !c write (iout,*) "NEWCORR BNEW2"
1609 !c write (iout,*) ((bnew2(ii,j,i),ii=1,3),j=1,2)
1611 read (ifourier,*,end=115,err=115) ccnew(kk,1,i)
1612 read (ifourier,*,end=115,err=115) ccnew(kk,2,i)
1614 !c write (iout,*) "NEWCORR CCNEW"
1615 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1617 read (ifourier,*,end=115,err=115) ddnew(kk,1,i)
1618 read (ifourier,*,end=115,err=115) ddnew(kk,2,i)
1620 !c write (iout,*) "NEWCORR DDNEW"
1621 !c write (iout,*) ((ddnew(ii,j,i),ii=1,3),j=1,2)
1625 read (ifourier,*,end=115,err=115) eenew(ii,jj,kk,i)
1629 !c write (iout,*) "NEWCORR EENEW1"
1630 !c write(iout,*)(((eenew(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1632 read (ifourier,*,end=115,err=115) e0new(ii,i)
1634 !c write (iout,*) (e0new(ii,i),ii=1,3)
1636 !c write (iout,*) "NEWCORR EENEW"
1639 ccnew(ii,1,i)=ccnew(ii,1,i)/2
1640 ccnew(ii,2,i)=ccnew(ii,2,i)/2
1641 ddnew(ii,1,i)=ddnew(ii,1,i)/2
1642 ddnew(ii,2,i)=ddnew(ii,2,i)/2
1647 bnew1(ii,1,-i)= bnew1(ii,1,i)
1648 bnew1(ii,2,-i)=-bnew1(ii,2,i)
1649 bnew2(ii,1,-i)= bnew2(ii,1,i)
1650 bnew2(ii,2,-i)=-bnew2(ii,2,i)
1653 !c ccnew(ii,1,i)=ccnew(ii,1,i)/2
1654 !c ccnew(ii,2,i)=ccnew(ii,2,i)/2
1655 !c ddnew(ii,1,i)=ddnew(ii,1,i)/2
1656 !c ddnew(ii,2,i)=ddnew(ii,2,i)/2
1657 ccnew(ii,1,-i)=ccnew(ii,1,i)
1658 ccnew(ii,2,-i)=-ccnew(ii,2,i)
1659 ddnew(ii,1,-i)=ddnew(ii,1,i)
1660 ddnew(ii,2,-i)=-ddnew(ii,2,i)
1662 e0new(1,-i)= e0new(1,i)
1663 e0new(2,-i)=-e0new(2,i)
1664 e0new(3,-i)=-e0new(3,i)
1666 eenew(kk,1,1,-i)= eenew(kk,1,1,i)
1667 eenew(kk,1,2,-i)=-eenew(kk,1,2,i)
1668 eenew(kk,2,1,-i)=-eenew(kk,2,1,i)
1669 eenew(kk,2,2,-i)= eenew(kk,2,2,i)
1673 write (iout,'(a)') "Coefficients of the multibody terms"
1674 do i=-nloctyp+1,nloctyp-1
1675 write (iout,*) "Type: ",onelet(iloctyp(i))
1676 write (iout,*) "Coefficients of the expansion of B1"
1678 write (iout,'(3hB1(,i1,1h),3f10.5)') j,(bnew1(k,j,i),k=1,3)
1680 write (iout,*) "Coefficients of the expansion of B2"
1682 write (iout,'(3hB2(,i1,1h),3f10.5)') j,(bnew2(k,j,i),k=1,3)
1684 write (iout,*) "Coefficients of the expansion of C"
1685 write (iout,'(3hC11,3f10.5)') (ccnew(j,1,i),j=1,3)
1686 write (iout,'(3hC12,3f10.5)') (ccnew(j,2,i),j=1,3)
1687 write (iout,*) "Coefficients of the expansion of D"
1688 write (iout,'(3hD11,3f10.5)') (ddnew(j,1,i),j=1,3)
1689 write (iout,'(3hD12,3f10.5)') (ddnew(j,2,i),j=1,3)
1690 write (iout,*) "Coefficients of the expansion of E"
1691 write (iout,'(2hE0,3f10.5)') (e0new(j,i),j=1,3)
1694 write (iout,'(1hE,2i1,2f10.5)') j,k,(eenew(l,j,k,i),l=1,2)
1699 IF (SPLIT_FOURIERTOR) THEN
1701 !c write (iout,*) "NEWCORR TOR",i
1702 read (ifourier,*,end=115,err=115)
1705 read (ifourier,*,end=115,err=115) bnew1tor(ii,j,i)
1708 !c write (iout,*) "NEWCORR BNEW1 TOR"
1709 !c write (iout,*) ((bnew1tor(ii,j,i),ii=1,3),j=1,2)
1712 read (ifourier,*,end=115,err=115) bnew2tor(ii,j,i)
1715 !c write (iout,*) "NEWCORR BNEW2 TOR"
1716 !c write (iout,*) ((bnew2tor(ii,j,i),ii=1,3),j=1,2)
1718 read (ifourier,*,end=115,err=115) ccnewtor(kk,1,i)
1719 read (ifourier,*,end=115,err=115) ccnewtor(kk,2,i)
1721 !c write (iout,*) "NEWCORR CCNEW TOR"
1722 !c write (iout,*) ((ccnew(ii,j,i),ii=1,3),j=1,2)
1724 read (ifourier,*,end=115,err=115) ddnewtor(kk,1,i)
1725 read (ifourier,*,end=115,err=115) ddnewtor(kk,2,i)
1727 !c write (iout,*) "NEWCORR DDNEW TOR"
1728 !c write (iout,*) ((ddnewtor(ii,j,i),ii=1,3),j=1,2)
1732 read (ifourier,*,end=115,err=115) eenewtor(ii,jj,kk,i)
1736 !c write (iout,*) "NEWCORR EENEW1 TOR"
1737 !c write(iout,*)(((eenewtor(ii,jj,kk,i),kk=1,2),jj=1,2),ii=1,2)
1739 read (ifourier,*,end=115,err=115) e0newtor(ii,i)
1741 !c write (iout,*) (e0newtor(ii,i),ii=1,3)
1743 !c write (iout,*) "NEWCORR EENEW TOR"
1746 ccnewtor(ii,1,i)=ccnewtor(ii,1,i)/2
1747 ccnewtor(ii,2,i)=ccnewtor(ii,2,i)/2
1748 ddnewtor(ii,1,i)=ddnewtor(ii,1,i)/2
1749 ddnewtor(ii,2,i)=ddnewtor(ii,2,i)/2
1754 bnew1tor(ii,1,-i)= bnew1tor(ii,1,i)
1755 bnew1tor(ii,2,-i)=-bnew1tor(ii,2,i)
1756 bnew2tor(ii,1,-i)= bnew2tor(ii,1,i)
1757 bnew2tor(ii,2,-i)=-bnew2tor(ii,2,i)
1760 ccnewtor(ii,1,-i)=ccnewtor(ii,1,i)
1761 ccnewtor(ii,2,-i)=-ccnewtor(ii,2,i)
1762 ddnewtor(ii,1,-i)=ddnewtor(ii,1,i)
1763 ddnewtor(ii,2,-i)=-ddnewtor(ii,2,i)
1765 e0newtor(1,-i)= e0newtor(1,i)
1766 e0newtor(2,-i)=-e0newtor(2,i)
1767 e0newtor(3,-i)=-e0newtor(3,i)
1769 eenewtor(kk,1,1,-i)= eenewtor(kk,1,1,i)
1770 eenewtor(kk,1,2,-i)=-eenewtor(kk,1,2,i)
1771 eenewtor(kk,2,1,-i)=-eenewtor(kk,2,1,i)
1772 eenewtor(kk,2,2,-i)= eenewtor(kk,2,2,i)
1776 write (iout,'(a)') &
1777 "Single-body coefficients of the torsional potentials"
1778 do i=-nloctyp+1,nloctyp-1
1779 write (iout,*) "Type: ",onelet(iloctyp(i))
1780 write (iout,*) "Coefficients of the expansion of B1tor"
1782 write (iout,'(3hB1(,i1,1h),3f10.5)') &
1783 j,(bnew1tor(k,j,i),k=1,3)
1785 write (iout,*) "Coefficients of the expansion of B2tor"
1787 write (iout,'(3hB2(,i1,1h),3f10.5)') &
1788 j,(bnew2tor(k,j,i),k=1,3)
1790 write (iout,*) "Coefficients of the expansion of Ctor"
1791 write (iout,'(3hC11,3f10.5)') (ccnewtor(j,1,i),j=1,3)
1792 write (iout,'(3hC12,3f10.5)') (ccnewtor(j,2,i),j=1,3)
1793 write (iout,*) "Coefficients of the expansion of Dtor"
1794 write (iout,'(3hD11,3f10.5)') (ddnewtor(j,1,i),j=1,3)
1795 write (iout,'(3hD12,3f10.5)') (ddnewtor(j,2,i),j=1,3)
1796 write (iout,*) "Coefficients of the expansion of Etor"
1797 write (iout,'(2hE0,3f10.5)') (e0newtor(j,i),j=1,3)
1800 write (iout,'(1hE,2i1,2f10.5)') &
1801 j,k,(eenewtor(l,j,k,i),l=1,2)
1807 do i=-nloctyp+1,nloctyp-1
1810 bnew1tor(ii,j,i)=bnew1(ii,j,i)
1815 bnew2tor(ii,j,i)=bnew2(ii,j,i)
1819 ccnewtor(ii,1,i)=ccnew(ii,1,i)
1820 ccnewtor(ii,2,i)=ccnew(ii,2,i)
1821 ddnewtor(ii,1,i)=ddnew(ii,1,i)
1822 ddnewtor(ii,2,i)=ddnew(ii,2,i)
1825 ENDIF !SPLIT_FOURIER_TOR
1827 allocate(ccold(2,2,-nloctyp-1:nloctyp+1))
1828 allocate(ddold(2,2,-nloctyp-1:nloctyp+1))
1829 allocate(eeold(2,2,-nloctyp-1:nloctyp+1))
1830 allocate(b(13,-nloctyp-1:nloctyp+1))
1832 write (iout,*) "Coefficients of the expansion of Eloc(l1,l2)"
1834 read (ifourier,*,end=115,err=115)
1835 read (ifourier,*,end=115,err=115) (b(ii,i),ii=1,13)
1837 write (iout,*) 'Type ',onelet(iloctyp(i))
1838 write (iout,'(a,i2,a,f10.5)') ('b(',ii,')=',b(ii,i),ii=1,13)
1852 !c B1tilde(1,i) = b(3)
1853 !c! B1tilde(2,i) =-b(5)
1854 !c! B1tilde(1,-i) =-b(3)
1855 !c! B1tilde(2,-i) =b(5)
1856 !c! b1tilde(1,i)=0.0d0
1857 !c b1tilde(2,i)=0.0d0
1862 !cc B1tilde(1,i) = b(3,i)
1863 !cc B1tilde(2,i) =-b(5,i)
1864 !c B1tilde(1,-i) =-b(3,i)
1865 !c B1tilde(2,-i) =b(5,i)
1866 !cc b1tilde(1,i)=0.0d0
1867 !cc b1tilde(2,i)=0.0d0
1868 !cc B2(1,i) = b(2,i)
1869 !cc B2(2,i) = b(4,i)
1871 !c B2(2,-i) =-b(4,i)
1875 CCold(1,1,i)= b(7,i)
1876 CCold(2,2,i)=-b(7,i)
1877 CCold(2,1,i)= b(9,i)
1878 CCold(1,2,i)= b(9,i)
1879 CCold(1,1,-i)= b(7,i)
1880 CCold(2,2,-i)=-b(7,i)
1881 CCold(2,1,-i)=-b(9,i)
1882 CCold(1,2,-i)=-b(9,i)
1887 !c Ctilde(1,1,i)= CCold(1,1,i)
1888 !c Ctilde(1,2,i)= CCold(1,2,i)
1889 !c Ctilde(2,1,i)=-CCold(2,1,i)
1890 !c Ctilde(2,2,i)=-CCold(2,2,i)
1895 !c Ctilde(1,1,i)= CCold(1,1,i)
1896 !c Ctilde(1,2,i)= CCold(1,2,i)
1897 !c Ctilde(2,1,i)=-CCold(2,1,i)
1898 !c Ctilde(2,2,i)=-CCold(2,2,i)
1900 !c Ctilde(1,1,i)=0.0d0
1901 !c Ctilde(1,2,i)=0.0d0
1902 !c Ctilde(2,1,i)=0.0d0
1903 !c Ctilde(2,2,i)=0.0d0
1904 DDold(1,1,i)= b(6,i)
1905 DDold(2,2,i)=-b(6,i)
1906 DDold(2,1,i)= b(8,i)
1907 DDold(1,2,i)= b(8,i)
1908 DDold(1,1,-i)= b(6,i)
1909 DDold(2,2,-i)=-b(6,i)
1910 DDold(2,1,-i)=-b(8,i)
1911 DDold(1,2,-i)=-b(8,i)
1916 !c Dtilde(1,1,i)= DD(1,1,i)
1917 !c Dtilde(1,2,i)= DD(1,2,i)
1918 !c Dtilde(2,1,i)=-DD(2,1,i)
1919 !c Dtilde(2,2,i)=-DD(2,2,i)
1921 !c Dtilde(1,1,i)=0.0d0
1922 !c Dtilde(1,2,i)=0.0d0
1923 !c Dtilde(2,1,i)=0.0d0
1924 !c Dtilde(2,2,i)=0.0d0
1925 EEold(1,1,i)= b(10,i)+b(11,i)
1926 EEold(2,2,i)=-b(10,i)+b(11,i)
1927 EEold(2,1,i)= b(12,i)-b(13,i)
1928 EEold(1,2,i)= b(12,i)+b(13,i)
1929 EEold(1,1,-i)= b(10,i)+b(11,i)
1930 EEold(2,2,-i)=-b(10,i)+b(11,i)
1931 EEold(2,1,-i)=-b(12,i)+b(13,i)
1932 EEold(1,2,-i)=-b(12,i)-b(13,i)
1933 write(iout,*) "TU DOCHODZE"
1939 !c ee(2,1,i)=ee(1,2,i)
1944 "Coefficients of the cumulants (independent of valence angles)"
1945 do i=-nloctyp+1,nloctyp-1
1946 write (iout,*) 'Type ',onelet(iloctyp(i))
1948 write(iout,'(2f10.5)') B(3,i),B(5,i)
1950 write(iout,'(2f10.5)') B(2,i),B(4,i)
1953 write (iout,'(2f10.5)') CCold(j,1,i),CCold(j,2,i)
1957 write (iout,'(2f10.5)') DDold(j,1,i),DDold(j,2,i)
1961 write (iout,'(2f10.5)') EEold(j,1,i),EEold(j,2,i)
1970 ! Read torsional parameters in old format
1972 allocate(itortyp(ntyp1)) !(-ntyp1:ntyp1)
1974 read (itorp,*,end=113,err=113) ntortyp,nterm_old
1975 if (lprint)write (iout,*) 'ntortyp,nterm',ntortyp,nterm_old
1976 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
1978 !el from energy module--------
1979 allocate(v1(nterm_old,ntortyp,ntortyp))
1980 allocate(v2(nterm_old,ntortyp,ntortyp)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor)
1981 !el---------------------------
1986 read (itorp,*,end=113,err=113) kk,v1(k,j,i),v2(k,j,i)
1992 write (iout,'(/a/)') 'Torsional constants:'
1995 write (iout,'(2i3,6f10.5)') i,j,(v1(k,i,j),k=1,nterm_old)
1996 write (iout,'(6x,6f10.5)') (v2(k,i,j),k=1,nterm_old)
2002 ! Read torsional parameters
2004 IF (TOR_MODE.eq.0) THEN
2005 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2006 read (itorp,*,end=113,err=113) ntortyp
2007 !el from energy module---------
2008 allocate(nterm(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2009 allocate(nlor(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2011 allocate(vlor1(maxlor,-ntortyp:ntortyp,-ntortyp:ntortyp)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2012 allocate(vlor2(maxlor,ntortyp,ntortyp))
2013 allocate(vlor3(maxlor,ntortyp,ntortyp)) !(maxlor,maxtor,maxtor)
2014 allocate(v0(-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2016 allocate(v1(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2))
2017 allocate(v2(maxterm,-ntortyp:ntortyp,-ntortyp:ntortyp,2)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2018 !el---------------------------
2021 !el---------------------------
2023 read (itorp,*,end=113,err=113) (itortyp(i),i=1,ntyp)
2025 itortyp(i)=-itortyp(-i)
2027 itortyp(ntyp1)=ntortyp
2028 itortyp(-ntyp1)=-ntortyp
2030 write (iout,*) 'ntortyp',ntortyp
2032 do j=-ntortyp+1,ntortyp-1
2033 read (itorp,*,end=113,err=113) nterm(i,j,iblock),&
2035 nterm(-i,-j,iblock)=nterm(i,j,iblock)
2036 nlor(-i,-j,iblock)=nlor(i,j,iblock)
2039 do k=1,nterm(i,j,iblock)
2040 read (itorp,*,end=113,err=113) kk,v1(k,i,j,iblock),&
2042 v1(k,-i,-j,iblock)=v1(k,i,j,iblock)
2043 v2(k,-i,-j,iblock)=-v2(k,i,j,iblock)
2044 v0ij=v0ij+si*v1(k,i,j,iblock)
2046 ! write(iout,*) i,j,k,iblock,nterm(i,j,iblock) !
2047 ! write(iout,*) v1(k,-i,-j,iblock),v1(k,i,j,iblock),&!
2048 ! v2(k,-i,-j,iblock),v2(k,i,j,iblock)!
2050 do k=1,nlor(i,j,iblock)
2051 read (itorp,*,end=113,err=113) kk,vlor1(k,i,j),&
2052 vlor2(k,i,j),vlor3(k,i,j)
2053 v0ij=v0ij+vlor1(k,i,j)/(1+vlor3(k,i,j)**2)
2056 v0(-i,-j,iblock)=v0ij
2062 write (iout,'(/a/)') 'Torsional constants:'
2064 do i=-ntortyp,ntortyp
2065 do j=-ntortyp,ntortyp
2066 write (iout,*) 'ityp',i,' jtyp',j
2067 write (iout,*) 'Fourier constants'
2068 do k=1,nterm(i,j,iblock)
2069 write (iout,'(2(1pe15.5))') v1(k,i,j,iblock),&
2072 write (iout,*) 'Lorenz constants'
2073 do k=1,nlor(i,j,iblock)
2074 write (iout,'(3(1pe15.5))') &
2075 vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2081 !elwrite (iout,'(/a/)') 'Torsional constants:',vlor1(k,i,j),vlor2(k,i,j),vlor3(k,i,j)
2083 ! 6/23/01 Read parameters for double torsionals
2085 !el from energy module------------
2086 allocate(v1c(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2087 allocate(v1s(2,maxtermd_1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2088 !(2,maxtermd_1,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2089 allocate(v2c(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2090 allocate(v2s(maxtermd_2,maxtermd_2,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2091 !(maxtermd_2,maxtermd_2,-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2092 allocate(ntermd_1(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2093 allocate(ntermd_2(-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,-ntortyp+1:ntortyp-1,2))
2094 !(-maxtor:maxtor,-maxtor:maxtor,-maxtor:maxtor,2)
2095 !---------------------------------
2099 do j=-ntortyp+1,ntortyp-1
2100 do k=-ntortyp+1,ntortyp-1
2101 read (itordp,'(3a1)',end=114,err=114) t1,t2,t3
2102 ! write (iout,*) "OK onelett",
2105 if (t1.ne.toronelet(i) .or. t2.ne.toronelet(j) &
2106 .or. t3.ne.toronelet(k)) then
2107 write (iout,*) "Error in double torsional parameter file",&
2110 call MPI_Finalize(Ierror)
2112 stop "Error in double torsional parameter file"
2114 read (itordp,*,end=114,err=114) ntermd_1(i,j,k,iblock),&
2115 ntermd_2(i,j,k,iblock)
2116 ntermd_1(-i,-j,-k,iblock)=ntermd_1(i,j,k,iblock)
2117 ntermd_2(-i,-j,-k,iblock)=ntermd_2(i,j,k,iblock)
2118 read (itordp,*,end=114,err=114) (v1c(1,l,i,j,k,iblock),l=1,&
2119 ntermd_1(i,j,k,iblock))
2120 read (itordp,*,end=114,err=114) (v1s(1,l,i,j,k,iblock),l=1,&
2121 ntermd_1(i,j,k,iblock))
2122 read (itordp,*,end=114,err=114) (v1c(2,l,i,j,k,iblock),l=1,&
2123 ntermd_1(i,j,k,iblock))
2124 read (itordp,*,end=114,err=114) (v1s(2,l,i,j,k,iblock),l=1,&
2125 ntermd_1(i,j,k,iblock))
2126 ! Martix of D parameters for one dimesional foureir series
2127 do l=1,ntermd_1(i,j,k,iblock)
2128 v1c(1,l,-i,-j,-k,iblock)=v1c(1,l,i,j,k,iblock)
2129 v1s(1,l,-i,-j,-k,iblock)=-v1s(1,l,i,j,k,iblock)
2130 v1c(2,l,-i,-j,-k,iblock)=v1c(2,l,i,j,k,iblock)
2131 v1s(2,l,-i,-j,-k,iblock)=-v1s(2,l,i,j,k,iblock)
2132 ! write(iout,*) "whcodze" ,
2133 ! & v1s(2,l,-i,-j,-k,iblock),v1s(2,l,i,j,k,iblock)
2135 read (itordp,*,end=114,err=114) ((v2c(l,m,i,j,k,iblock),&
2136 v2c(m,l,i,j,k,iblock),v2s(l,m,i,j,k,iblock),&
2137 v2s(m,l,i,j,k,iblock),&
2138 m=1,l-1),l=1,ntermd_2(i,j,k,iblock))
2139 ! Martix of D parameters for two dimesional fourier series
2140 do l=1,ntermd_2(i,j,k,iblock)
2142 v2c(l,m,-i,-j,-k,iblock)=v2c(l,m,i,j,k,iblock)
2143 v2c(m,l,-i,-j,-k,iblock)=v2c(m,l,i,j,k,iblock)
2144 v2s(l,m,-i,-j,-k,iblock)=-v2s(l,m,i,j,k,iblock)
2145 v2s(m,l,-i,-j,-k,iblock)=-v2s(m,l,i,j,k,iblock)
2154 write (iout,*) 'Constants for double torsionals'
2157 do j=-ntortyp+1,ntortyp-1
2158 do k=-ntortyp+1,ntortyp-1
2159 write (iout,*) 'ityp',i,' jtyp',j,' ktyp',k,&
2160 ' nsingle',ntermd_1(i,j,k,iblock),&
2161 ' ndouble',ntermd_2(i,j,k,iblock)
2163 write (iout,*) 'Single angles:'
2164 do l=1,ntermd_1(i,j,k,iblock)
2165 write (iout,'(i5,2f10.5,5x,2f10.5,5x,2f10.5)') l,&
2166 v1c(1,l,i,j,k,iblock),v1s(1,l,i,j,k,iblock),&
2167 v1c(2,l,i,j,k,iblock),v1s(2,l,i,j,k,iblock),&
2168 v1s(1,l,-i,-j,-k,iblock),v1s(2,l,-i,-j,-k,iblock)
2171 write (iout,*) 'Pairs of angles:'
2172 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2173 do l=1,ntermd_2(i,j,k,iblock)
2174 write (iout,'(i5,20f10.5)') &
2175 l,(v2c(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock))
2178 write (iout,'(3x,20i10)') (l,l=1,ntermd_2(i,j,k,iblock))
2179 do l=1,ntermd_2(i,j,k,iblock)
2180 write (iout,'(i5,20f10.5)') &
2181 l,(v2s(l,m,i,j,k,iblock),m=1,ntermd_2(i,j,k,iblock)),&
2182 (v2s(l,m,-i,-j,-k,iblock),m=1,ntermd_2(i,j,k,iblock))
2190 ELSE IF (TOR_MODE.eq.1) THEN
2192 !C read valence-torsional parameters
2193 read (itorp,*,end=121,err=121) ntortyp
2195 write (iout,*) "Valence-torsional parameters read in ntortyp",&
2197 read (itorp,*,end=121,err=121) (itortyp(i),i=1,ntyp)
2198 write (iout,*) "itortyp_kcc",(itortyp(i),i=1,ntyp)
2201 itype2loc(i)=itortyp(i)
2205 itortyp(i)=-itortyp(-i)
2207 do i=-ntortyp+1,ntortyp-1
2208 do j=-ntortyp+1,ntortyp-1
2209 !C first we read the cos and sin gamma parameters
2210 read (itorp,'(13x,a)',end=121,err=121) string
2211 write (iout,*) i,j,string
2212 read (itorp,*,end=121,err=121) &
2213 nterm_kcc(j,i),nterm_kcc_Tb(j,i)
2214 !C read (itorkcc,*,end=121,err=121) nterm_kcc_Tb(j,i)
2215 do k=1,nterm_kcc(j,i)
2216 do l=1,nterm_kcc_Tb(j,i)
2217 do ll=1,nterm_kcc_Tb(j,i)
2218 read (itorp,*,end=121,err=121) ii,jj,kk, &
2219 v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2227 !c AL 4/8/16: Calculate coefficients from one-body parameters
2229 allocate(itortyp(-ntyp1:ntyp1)) !(-ntyp1:ntyp1)
2230 allocate(nterm_kcc(-ntyp1:ntyp1,-ntyp1:ntyp1))
2231 allocate(nterm_kcc_Tb(-ntyp1:ntyp1,-ntyp1:ntyp1))
2232 allocate(v1_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2233 allocate(v2_kcc(6,6,6,-ntyp1:ntyp1,-ntyp1:ntyp1))
2236 print *,i,itortyp(i)
2237 itortyp(i)=itype2loc(i)
2240 "Val-tor parameters calculated from cumulant coefficients ntortyp"&
2242 do i=-ntortyp+1,ntortyp-1
2243 do j=-ntortyp+1,ntortyp-1
2246 do k=1,nterm_kcc_Tb(j,i)
2247 do l=1,nterm_kcc_Tb(j,i)
2248 v1_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,1,j)&
2249 +bnew1tor(k,2,i)*bnew2tor(l,2,j)
2250 v2_kcc(k,l,1,i,j)=bnew1tor(k,1,i)*bnew2tor(l,2,j)&
2251 +bnew1tor(k,2,i)*bnew2tor(l,1,j)
2254 do k=1,nterm_kcc_Tb(j,i)
2255 do l=1,nterm_kcc_Tb(j,i)
2257 v1_kcc(k,l,2,i,j)=-(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2258 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2259 v2_kcc(k,l,2,i,j)=-(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2260 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2262 v1_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,1,i)*ddnewtor(l,1,j) &
2263 -ccnewtor(k,2,i)*ddnewtor(l,2,j))
2264 v2_kcc(k,l,2,i,j)=-0.25*(ccnewtor(k,2,i)*ddnewtor(l,1,j) &
2265 +ccnewtor(k,1,i)*ddnewtor(l,2,j))
2269 !c f(theta,gamma)=-(b21(theta)*b11(theta)+b12(theta)*b22(theta))*cos(gamma)-(b22(theta)*b11(theta)+b21(theta)*b12(theta))*sin(gamma)+(c11(theta)*d11(theta)-c12(theta)*d12(theta))*cos(2*gamma)+(c12(theta)*d11(theta)+c11(theta)*d12(theta))*sin(2*gamma)
2273 write (iout,*) "TOR_MODE>1 only with NEWCORR"
2278 if (tor_mode.gt.0 .and. lprint) then
2279 !c Print valence-torsional parameters
2280 write (iout,'(a)') &
2281 "Parameters of the valence-torsional potentials"
2282 do i=-ntortyp+1,ntortyp-1
2283 do j=-ntortyp+1,ntortyp-1
2284 write (iout,'(3a)') "Type ",toronelet(i),toronelet(j)
2285 write (iout,'(3a5,2a15)') "itor","ival","jval","v_kcc","v2_kcc"
2286 do k=1,nterm_kcc(j,i)
2287 do l=1,nterm_kcc_Tb(j,i)
2288 do ll=1,nterm_kcc_Tb(j,i)
2289 write (iout,'(3i5,2f15.4)')&
2290 k,l-1,ll-1,v1_kcc(ll,l,k,j,i),v2_kcc(ll,l,k,j,i)
2298 allocate(itortyp_nucl(ntyp1_molec(2))) !(-ntyp1:ntyp1)
2299 read (itorp_nucl,*,end=113,err=113) ntortyp_nucl
2300 ! print *,"ntortyp_nucl",ntortyp_nucl,ntyp_molec(2)
2301 !el from energy module---------
2302 allocate(nterm_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2303 allocate(nlor_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2305 allocate(vlor1_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,-maxtor:maxtor,-maxtor:maxtor)
2306 allocate(vlor2_nucl(maxlor,ntortyp_nucl,ntortyp_nucl))
2307 allocate(vlor3_nucl(maxlor,ntortyp_nucl,ntortyp_nucl)) !(maxlor,maxtor,maxtor)
2308 allocate(v0_nucl(ntortyp_nucl,ntortyp_nucl)) !(-maxtor:maxtor,-maxtor:maxtor,2)
2310 allocate(v1_nucl(maxterm,ntortyp_nucl,ntortyp_nucl))
2311 allocate(v2_nucl(maxterm,ntortyp_nucl,ntortyp_nucl)) !(maxterm,-maxtor:maxtor,-maxtor:maxtor,2)
2312 !el---------------------------
2315 !el--------------------
2316 read (itorp_nucl,*,end=113,err=113) &
2317 (itortyp_nucl(i),i=1,ntyp_molec(2))
2318 ! print *,itortyp_nucl(:)
2319 !c write (iout,*) 'ntortyp',ntortyp
2322 read (itorp_nucl,*,end=113,err=113) nterm_nucl(i,j),nlor_nucl(i,j)
2323 ! print *,nterm_nucl(i,j),nlor_nucl(i,j)
2326 do k=1,nterm_nucl(i,j)
2327 read (itorp_nucl,*,end=113,err=113) kk,v1_nucl(k,i,j),v2_nucl(k,i,j)
2328 v0ij=v0ij+si*v1_nucl(k,i,j)
2331 do k=1,nlor_nucl(i,j)
2332 read (itorp,*,end=113,err=113) kk,vlor1_nucl(k,i,j),&
2333 vlor2_nucl(k,i,j),vlor3_nucl(k,i,j)
2334 v0ij=v0ij+vlor1_nucl(k,i,j)/(1+vlor3_nucl(k,i,j)**2)
2340 ! Read of Side-chain backbone correlation parameters
2341 ! Modified 11 May 2012 by Adasko
2344 read (isccor,*,end=119,err=119) nsccortyp
2346 !el from module energy-------------
2347 allocate(nlor_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2348 allocate(vlor1sccor(maxterm_sccor,nsccortyp,nsccortyp))
2349 allocate(vlor2sccor(maxterm_sccor,nsccortyp,nsccortyp))
2350 allocate(vlor3sccor(maxterm_sccor,nsccortyp,nsccortyp)) !(maxterm_sccor,20,20)
2351 !-----------------------------------
2353 !el from module energy-------------
2354 allocate(isccortyp(-ntyp:ntyp)) !(-ntyp:ntyp)
2356 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2358 isccortyp(i)=-isccortyp(-i)
2360 iscprol=isccortyp(20)
2361 ! write (iout,*) 'ntortyp',ntortyp
2363 !c maxinter is maximum interaction sites
2364 !el from module energy---------
2365 allocate(nterm_sccor(-nsccortyp:nsccortyp,-nsccortyp:nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2366 allocate(v1sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2367 -nsccortyp:nsccortyp))
2368 allocate(v2sccor(maxterm_sccor,maxinter,-nsccortyp:nsccortyp,&
2369 -nsccortyp:nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2370 allocate(v0sccor(maxinter,-nsccortyp:nsccortyp,&
2371 -nsccortyp:nsccortyp)) !(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2372 !-----------------------------------
2376 read (isccor,*,end=119,err=119) &
2377 nterm_sccor(i,j),nlor_sccor(i,j)
2383 nterm_sccor(-i,j)=nterm_sccor(i,j)
2384 nterm_sccor(-i,-j)=nterm_sccor(i,j)
2385 nterm_sccor(i,-j)=nterm_sccor(i,j)
2386 do k=1,nterm_sccor(i,j)
2387 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2389 if (j.eq.iscprol) then
2390 if (i.eq.isccortyp(10)) then
2391 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2392 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2394 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)*0.5d0 &
2395 +v2sccor(k,l,i,j)*dsqrt(0.75d0)
2396 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)*0.5d0 &
2397 +v1sccor(k,l,i,j)*dsqrt(0.75d0)
2398 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2399 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2400 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2401 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2404 if (i.eq.isccortyp(10)) then
2405 v1sccor(k,l,i,-j)=v1sccor(k,l,i,j)
2406 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2408 if (j.eq.isccortyp(10)) then
2409 v1sccor(k,l,-i,j)=v1sccor(k,l,i,j)
2410 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,j)
2412 v1sccor(k,l,i,-j)=-v1sccor(k,l,i,j)
2413 v2sccor(k,l,i,-j)=-v2sccor(k,l,i,j)
2414 v1sccor(k,l,-i,-j)=v1sccor(k,l,i,j)
2415 v2sccor(k,l,-i,-j)=-v2sccor(k,l,i,j)
2416 v1sccor(k,l,-i,j)=v1sccor(k,l,i,-j)
2417 v2sccor(k,l,-i,j)=-v2sccor(k,l,i,-j)
2421 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2422 v0ijsccor1=v0ijsccor+si*v1sccor(k,l,-i,j)
2423 v0ijsccor2=v0ijsccor+si*v1sccor(k,l,i,-j)
2424 v0ijsccor3=v0ijsccor+si*v1sccor(k,l,-i,-j)
2427 do k=1,nlor_sccor(i,j)
2428 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2429 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2430 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2431 (1+vlor3sccor(k,i,j)**2)
2433 v0sccor(l,i,j)=v0ijsccor
2434 v0sccor(l,-i,j)=v0ijsccor1
2435 v0sccor(l,i,-j)=v0ijsccor2
2436 v0sccor(l,-i,-j)=v0ijsccor3
2442 !el from module energy-------------
2443 allocate(isccortyp(ntyp)) !(-ntyp:ntyp)
2445 read (isccor,*,end=119,err=119) (isccortyp(i),i=1,ntyp)
2446 ! write (iout,*) 'ntortyp',ntortyp
2448 !c maxinter is maximum interaction sites
2449 !el from module energy---------
2450 allocate(nterm_sccor(nsccortyp,nsccortyp)) !(-ntyp:ntyp,-ntyp:ntyp)
2451 allocate(v1sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp))
2452 allocate(v2sccor(maxterm_sccor,maxinter,nsccortyp,nsccortyp)) !(maxterm_sccor,3,-ntyp:ntyp,-ntyp:ntyp)
2453 allocate(v0sccor(maxinter,nsccortyp,nsccortyp)) !???(maxterm_sccor,-ntyp:ntyp,-ntyp:ntyp)
2454 !-----------------------------------
2458 read (isccor,*,end=119,err=119) &
2459 nterm_sccor(i,j),nlor_sccor(i,j)
2463 do k=1,nterm_sccor(i,j)
2464 read (isccor,*,end=119,err=119) kk,v1sccor(k,l,i,j),&
2466 v0ijsccor=v0ijsccor+si*v1sccor(k,l,i,j)
2469 do k=1,nlor_sccor(i,j)
2470 read (isccor,*,end=119,err=119) kk,vlor1sccor(k,i,j),&
2471 vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2472 v0ijsccor=v0ijsccor+vlor1sccor(k,i,j)/ &
2473 (1+vlor3sccor(k,i,j)**2)
2475 v0sccor(l,i,j)=v0ijsccor !el ,iblock
2484 write (iout,'(/a/)') 'Torsional constants:'
2487 write (iout,*) 'ityp',i,' jtyp',j
2488 write (iout,*) 'Fourier constants'
2489 do k=1,nterm_sccor(i,j)
2490 write (iout,'(2(1pe15.5))') v1sccor(k,l,i,j),v2sccor(k,l,i,j)
2492 write (iout,*) 'Lorenz constants'
2493 do k=1,nlor_sccor(i,j)
2494 write (iout,'(3(1pe15.5))') &
2495 vlor1sccor(k,i,j),vlor2sccor(k,i,j),vlor3sccor(k,i,j)
2502 ! 9/18/99 (AL) Read coefficients of the Fourier expansion of the local
2503 ! interaction energy of the Gly, Ala, and Pro prototypes.
2506 ! Read electrostatic-interaction parameters
2511 write (iout,'(/a)') 'Electrostatic interaction constants:'
2512 write (iout,'(1x,a,1x,a,10x,a,11x,a,11x,a,11x,a)') &
2513 'IT','JT','APP','BPP','AEL6','AEL3'
2515 read (ielep,*,end=116,err=116) ((epp(i,j),j=1,2),i=1,2)
2516 read (ielep,*,end=116,err=116) ((rpp(i,j),j=1,2),i=1,2)
2517 read (ielep,*,end=116,err=116) ((elpp6(i,j),j=1,2),i=1,2)
2518 read (ielep,*,end=116,err=116) ((elpp3(i,j),j=1,2),i=1,2)
2523 app (i,j)=epp(i,j)*rri*rri
2524 bpp (i,j)=-2.0D0*epp(i,j)*rri
2525 ael6(i,j)=elpp6(i,j)*4.2D0**6
2526 ael3(i,j)=elpp3(i,j)*4.2D0**3
2528 if (lprint) write(iout,'(2i3,4(1pe15.4))')i,j,app(i,j),bpp(i,j),&
2534 ! Read side-chain interaction parameters.
2536 !el from module energy - COMMON.INTERACT-------
2537 allocate(eps(ntyp,ntyp),sigmaii(ntyp,ntyp),rs0(ntyp,ntyp)) !(ntyp,ntyp)
2538 allocate(augm(ntyp,ntyp)) !(ntyp,ntyp)
2539 allocate(eps_scp(ntyp,2),rscp(ntyp,2)) !(ntyp,2)
2541 allocate(sigma0(ntyp),rr0(ntyp),sigii(ntyp)) !(ntyp)
2542 allocate(chip(ntyp1),alp(ntyp1)) !(ntyp)
2543 allocate(epslip(ntyp,ntyp))
2551 !--------------------------------
2553 read (isidep,*,end=117,err=117) ipot,expon
2554 if (ipot.lt.1 .or. ipot.gt.5) then
2555 write (iout,'(2a)') 'Error while reading SC interaction',&
2556 'potential file - unknown potential type.'
2558 call MPI_Finalize(Ierror)
2564 write(iout,'(/3a,2i3)') 'Potential is ',potname(ipot),&
2565 ', exponents are ',expon,2*expon
2566 ! goto (10,20,30,30,40) ipot
2568 !----------------------- LJ potential ---------------------------------
2570 ! 10 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2571 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2572 (sigma0(i),i=1,ntyp)
2574 write (iout,'(/a/)') 'Parameters of the LJ potential:'
2575 write (iout,'(a/)') 'The epsilon array:'
2576 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2577 write (iout,'(/a)') 'One-body parameters:'
2578 write (iout,'(a,4x,a)') 'residue','sigma'
2579 write (iout,'(a3,6x,f10.5)') (restyp(i,1),sigma0(i),i=1,ntyp)
2582 !----------------------- LJK potential --------------------------------
2584 ! 20 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2585 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2586 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp)
2588 write (iout,'(/a/)') 'Parameters of the LJK potential:'
2589 write (iout,'(a/)') 'The epsilon array:'
2590 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2591 write (iout,'(/a)') 'One-body parameters:'
2592 write (iout,'(a,4x,2a)') 'residue',' sigma ',' r0 '
2593 write (iout,'(a3,6x,2f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2597 !---------------------- GB or BP potential -----------------------------
2600 ! print *,"I AM in SCELE",scelemode
2601 if (scelemode.eq.0) then
2603 read (isidep,*,end=117,err=117)(eps(i,j),j=i,ntyp)
2605 read (isidep,*,end=117,err=117)(sigma0(i),i=1,ntyp)
2606 read (isidep,*,end=117,err=117)(sigii(i),i=1,ntyp)
2607 read (isidep,*,end=117,err=117)(chip(i),i=1,ntyp)
2608 read (isidep,*,end=117,err=117)(alp(i),i=1,ntyp)
2610 read (isidep,*,end=117,err=117)(epslip(i,j),j=i,ntyp)
2613 ! For the GB potential convert sigma'**2 into chi'
2616 chip(i)=(chip(i)-1.0D0)/(chip(i)+1.0D0)
2620 write (iout,'(/a/)') 'Parameters of the BP potential:'
2621 write (iout,'(a/)') 'The epsilon array:'
2622 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2623 write (iout,'(/a)') 'One-body parameters:'
2624 write (iout,'(a,4x,4a)') 'residue',' sigma ','s||/s_|_^2',&
2626 write (iout,'(a3,6x,4f10.5)') (restyp(i,1),sigma0(i),sigii(i),&
2627 chip(i),alp(i),i=1,ntyp)
2630 ! print *,ntyp,"NTYP"
2631 allocate(icharge(ntyp1))
2632 ! print *,ntyp,icharge(i)
2634 read (isidep,*) (icharge(i),i=1,ntyp)
2635 print *,ntyp,icharge(i)
2636 ! if(.not.allocated(eps)) allocate(eps(-ntyp
2637 !c write (2,*) "icharge",(icharge(i),i=1,ntyp)
2638 allocate(alphapol(ntyp,ntyp),epshead(ntyp,ntyp),sig0head(ntyp,ntyp))
2639 allocate(sigiso1(ntyp,ntyp),rborn(ntyp,ntyp),sigmap1(ntyp,ntyp))
2640 allocate(sigmap2(ntyp,ntyp),sigiso2(ntyp,ntyp))
2641 allocate(chis(ntyp,ntyp),wquad(ntyp,ntyp),chipp(ntyp,ntyp))
2642 allocate(epsintab(ntyp,ntyp))
2643 allocate(dtail(2,ntyp,ntyp))
2644 allocate(alphasur(4,ntyp,ntyp),alphiso(4,ntyp,ntyp))
2645 allocate(wqdip(2,ntyp,ntyp))
2646 allocate(wstate(4,ntyp,ntyp))
2647 allocate(dhead(2,2,ntyp,ntyp))
2648 allocate(nstate(ntyp,ntyp))
2649 allocate(debaykap(ntyp,ntyp))
2651 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2652 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2656 ! write (*,*) "Im in ALAB", i, " ", j
2658 eps(i,j),sigma(i,j),chi(i,j),chi(j,i),chipp(i,j),chipp(j,i), &
2659 (alphasur(k,i,j),k=1,4),sigmap1(i,j),sigmap2(i,j), &
2660 chis(i,j),chis(j,i), &
2661 nstate(i,j),(wstate(k,i,j),k=1,4), &
2662 dhead(1,1,i,j),dhead(1,2,i,j),dhead(2,1,i,j),dhead(2,2,i,j),&
2663 dtail(1,i,j),dtail(2,i,j), &
2664 epshead(i,j),sig0head(i,j), &
2665 rborn(i,j),rborn(j,i),(wqdip(k,i,j),k=1,2),wquad(i,j), &
2666 alphapol(i,j),alphapol(j,i), &
2667 (alphiso(k,i,j),k=1,4),sigiso1(i,j),sigiso2(i,j),epsintab(i,j),debaykap(i,j)
2668 ! print *,eps(i,j),sigma(i,j),"SIGMAP",i,j,sigmap1(i,j),sigmap2(j,i)
2674 sigma(i,j) = sigma(j,i)
2675 sigmap1(i,j)=sigmap1(j,i)
2676 sigmap2(i,j)=sigmap2(j,i)
2677 sigiso1(i,j)=sigiso1(j,i)
2678 sigiso2(i,j)=sigiso2(j,i)
2679 ! print *,"ATU",sigma(j,i),sigma(i,j),i,j
2680 nstate(i,j) = nstate(j,i)
2681 dtail(1,i,j) = dtail(1,j,i)
2682 dtail(2,i,j) = dtail(2,j,i)
2684 alphasur(k,i,j) = alphasur(k,j,i)
2685 wstate(k,i,j) = wstate(k,j,i)
2686 alphiso(k,i,j) = alphiso(k,j,i)
2689 dhead(2,1,i,j) = dhead(1,1,j,i)
2690 dhead(2,2,i,j) = dhead(1,2,j,i)
2691 dhead(1,1,i,j) = dhead(2,1,j,i)
2692 dhead(1,2,i,j) = dhead(2,2,j,i)
2694 epshead(i,j) = epshead(j,i)
2695 sig0head(i,j) = sig0head(j,i)
2698 wqdip(k,i,j) = wqdip(k,j,i)
2701 wquad(i,j) = wquad(j,i)
2702 epsintab(i,j) = epsintab(j,i)
2703 debaykap(i,j)=debaykap(j,i)
2704 ! if (epsintab(i,j).ne.1.0) print *,"WHAT?",i,j,epsintab(i,j)
2709 !--------------------- GBV potential -----------------------------------
2711 ! 40 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2712 read (isidep,*,end=117,err=117)((eps(i,j),j=i,ntyp),i=1,ntyp),&
2713 (sigma0(i),i=1,ntyp),(rr0(i),i=1,ntyp),(sigii(i),i=1,ntyp),&
2714 (chip(i),i=1,ntyp),(alp(i),i=1,ntyp)
2716 write (iout,'(/a/)') 'Parameters of the GBV potential:'
2717 write (iout,'(a/)') 'The epsilon array:'
2718 call printmat(ntyp,ntyp,ntyp,iout,restyp,eps)
2719 write (iout,'(/a)') 'One-body parameters:'
2720 write (iout,'(a,4x,5a)') 'residue',' sigma ',' r0 ',&
2721 's||/s_|_^2',' chip ',' alph '
2722 write (iout,'(a3,6x,5f10.5)') (restyp(i,1),sigma0(i),rr0(i),&
2723 sigii(i),chip(i),alp(i),i=1,ntyp)
2726 write(iout,*)"Wrong ipot"
2732 !-----------------------------------------------------------------------
2733 ! Calculate the "working" parameters of SC interactions.
2735 !el from module energy - COMMON.INTERACT-------
2736 allocate(aa_aq(ntyp1,ntyp1),bb_aq(ntyp1,ntyp1))
2737 if (.not.allocated(chi)) allocate(chi(ntyp1,ntyp1)) !(ntyp,ntyp)
2738 allocate(aa_lip(ntyp1,ntyp1),bb_lip(ntyp1,ntyp1)) !(ntyp,ntyp)
2739 if (.not.allocated(sigma)) allocate(sigma(0:ntyp1,0:ntyp1))
2740 allocate(r0(ntyp1,ntyp1)) !(0:ntyp1,0:ntyp1)
2741 allocate(acavtub(ntyp1),bcavtub(ntyp1),ccavtub(ntyp1),&
2743 allocate(sc_aa_tube_par(ntyp1),sc_bb_tube_par(ntyp1),&
2749 if (scelemode.eq.0) then
2758 sc_aa_tube_par(:)=0.0d0
2759 sc_bb_tube_par(:)=0.0d0
2761 !--------------------------------
2766 epslip(i,j)=epslip(j,i)
2769 if (scelemode.eq.0) then
2772 sigma(i,j)=dsqrt(sigma0(i)**2+sigma0(j)**2)
2773 sigma(j,i)=sigma(i,j)
2774 rs0(i,j)=dwa16*sigma(i,j)
2779 if (lprint) write (iout,'(/a/10x,7a/72(1h-))') &
2780 'Working parameters of the SC interactions:',&
2781 ' a ',' b ',' augm ',' sigma ',' r0 ',&
2786 if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2788 ! print *,"SIGMA ZLA?",sigma(i,j)
2796 sigeps=dsign(1.0D0,epsij)
2798 aa_aq(i,j)=epsij*rrij*rrij
2799 ! print *,"ADASKO",epsij,rrij,expon
2800 bb_aq(i,j)=-sigeps*epsij*rrij
2801 aa_aq(j,i)=aa_aq(i,j)
2802 bb_aq(j,i)=bb_aq(i,j)
2803 epsijlip=epslip(i,j)
2804 sigeps=dsign(1.0D0,epsijlip)
2805 epsijlip=dabs(epsijlip)
2806 aa_lip(i,j)=epsijlip*rrij*rrij
2807 bb_lip(i,j)=-sigeps*epsijlip*rrij
2808 aa_lip(j,i)=aa_lip(i,j)
2809 bb_lip(j,i)=bb_lip(i,j)
2810 !C write(iout,*) aa_lip
2811 if ((ipot.gt.2).and. (scelemode.eq.0)) then
2812 sigt1sq=sigma0(i)**2
2813 sigt2sq=sigma0(j)**2
2816 ratsig1=sigt2sq/sigt1sq
2817 ratsig2=1.0D0/ratsig1
2818 chi(i,j)=(sigii1-1.0D0)/(sigii1+ratsig1)
2819 if (j.gt.i) chi(j,i)=(sigii2-1.0D0)/(sigii2+ratsig2)
2820 rsum_max=dsqrt(sigii1*sigt1sq+sigii2*sigt2sq)
2824 ! if (ipot.eq.1 .or. ipot.eq.3 .or. ipot.eq.4) then
2825 sigmaii(i,j)=rsum_max
2826 sigmaii(j,i)=rsum_max
2828 ! sigmaii(i,j)=r0(i,j)
2829 ! sigmaii(j,i)=r0(i,j)
2831 !d write (iout,*) i,j,r0(i,j),sigma(i,j),rsum_max
2832 if ((ipot.eq.2 .or. ipot.eq.5) .and. r0(i,j).gt.rsum_max) then
2833 r_augm=sigma(i,j)*(rrij-sigma(i,j))/rrij
2834 augm(i,j)=epsij*r_augm**(2*expon)
2835 ! augm(i,j)=0.5D0**(2*expon)*aa(i,j)
2842 write (iout,'(2(a3,2x),3(1pe10.3),5(0pf8.3))') &
2843 restyp(i,1),restyp(j,1),aa_aq(i,j),bb_aq(i,j),augm(i,j),&
2844 sigma(i,j),r0(i,j),chi(i,j),chi(j,i)
2849 allocate(eps_nucl(ntyp_molec(2),ntyp_molec(2)))
2850 allocate(sigma_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2851 allocate(elpp6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2852 allocate(elpp3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2853 allocate(elpp63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2854 allocate(elpp32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2855 allocate(chi_nucl(ntyp_molec(2),ntyp_molec(2)),chip_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2856 allocate(ael3_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2857 allocate(ael6_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2858 allocate(ael32_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2859 allocate(ael63_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2860 allocate(aa_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2861 allocate(bb_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2862 allocate(r0_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp,2)
2863 allocate(sigmaii_nucl(ntyp_molec(2),ntyp_molec(2))) !(ntyp_molec(2),ntyp_molec(2))
2864 allocate(eps_scp_nucl(ntyp_molec(2)),rscp_nucl(ntyp_molec(2))) !(ntyp,2)
2873 read (isidep_nucl,*) ipot_nucl
2874 ! print *,"TU?!",ipot_nucl
2875 if (ipot_nucl.eq.1) then
2876 do i=1,ntyp_molec(2)
2877 do j=i,ntyp_molec(2)
2878 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),elpp6_nucl(i,j),&
2879 elpp3_nucl(i,j), elpp63_nucl(i,j),elpp32_nucl(i,j)
2883 do i=1,ntyp_molec(2)
2884 do j=i,ntyp_molec(2)
2885 read (isidep_nucl,*) eps_nucl(i,j),sigma_nucl(i,j),chi_nucl(i,j),&
2886 chi_nucl(j,i),chip_nucl(i,j),chip_nucl(j,i),&
2887 elpp6_nucl(i,j),elpp3_nucl(i,j),elpp63_nucl(i,j),elpp32_nucl(i,j)
2891 ! rpp(1,1)=2**(1.0/6.0)*5.16158
2892 do i=1,ntyp_molec(2)
2893 do j=i,ntyp_molec(2)
2894 rrij=sigma_nucl(i,j)
2898 epsij=4*eps_nucl(i,j)
2899 sigeps=dsign(1.0D0,epsij)
2901 aa_nucl(i,j)=epsij*rrij*rrij
2902 bb_nucl(i,j)=-sigeps*epsij*rrij
2903 ael3_nucl(i,j)=elpp3_nucl(i,j)*dsqrt(rrij)
2904 ael6_nucl(i,j)=elpp6_nucl(i,j)*rrij
2905 ael63_nucl(i,j)=elpp63_nucl(i,j)*rrij
2906 ael32_nucl(i,j)=elpp32_nucl(i,j)*rrij
2907 sigmaii_nucl(i,j)=sigma_nucl(i,j)/sqrt(1-(chi_nucl(i,j)+chi_nucl(j,i)- &
2908 2*chi_nucl(i,j)*chi_nucl(j,i))/(1-chi_nucl(i,j)*chi_nucl(j,i)))
2911 aa_nucl(i,j)=aa_nucl(j,i)
2912 bb_nucl(i,j)=bb_nucl(j,i)
2913 ael3_nucl(i,j)=ael3_nucl(j,i)
2914 ael6_nucl(i,j)=ael6_nucl(j,i)
2915 ael63_nucl(i,j)=ael63_nucl(j,i)
2916 ael32_nucl(i,j)=ael32_nucl(j,i)
2917 elpp3_nucl(i,j)=elpp3_nucl(j,i)
2918 elpp6_nucl(i,j)=elpp6_nucl(j,i)
2919 elpp63_nucl(i,j)=elpp63_nucl(j,i)
2920 elpp32_nucl(i,j)=elpp32_nucl(j,i)
2921 eps_nucl(i,j)=eps_nucl(j,i)
2922 sigma_nucl(i,j)=sigma_nucl(j,i)
2923 sigmaii_nucl(i,j)=sigmaii_nucl(j,i)
2927 write(iout,*) "tube param"
2928 read(itube,*) epspeptube,sigmapeptube,acavtubpep,bcavtubpep, &
2929 ccavtubpep,dcavtubpep,tubetranenepep
2930 sigmapeptube=sigmapeptube**6
2931 sigeps=dsign(1.0D0,epspeptube)
2932 epspeptube=dabs(epspeptube)
2933 pep_aa_tube=4.0d0*epspeptube*sigmapeptube**2
2934 pep_bb_tube=-sigeps*4.0d0*epspeptube*sigmapeptube
2935 write(iout,*) pep_aa_tube,pep_bb_tube,tubetranenepep
2937 read(itube,*) epssctube,sigmasctube,acavtub(i),bcavtub(i), &
2938 ccavtub(i),dcavtub(i),tubetranene(i)
2939 sigmasctube=sigmasctube**6
2940 sigeps=dsign(1.0D0,epssctube)
2941 epssctube=dabs(epssctube)
2942 sc_aa_tube_par(i)=4.0d0*epssctube*sigmasctube**2
2943 sc_bb_tube_par(i)=-sigeps*4.0d0*epssctube*sigmasctube
2944 write(iout,*) sc_aa_tube_par(i), sc_bb_tube_par(i),tubetranene(i)
2946 !-----------------READING SC BASE POTENTIALS-----------------------------
2947 allocate(eps_scbase(ntyp_molec(1),ntyp_molec(2)))
2948 allocate(sigma_scbase(ntyp_molec(1),ntyp_molec(2)))
2949 allocate(chi_scbase(ntyp_molec(1),ntyp_molec(2),2))
2950 allocate(chipp_scbase(ntyp_molec(1),ntyp_molec(2),2))
2951 allocate(alphasur_scbase(4,ntyp_molec(1),ntyp_molec(2)))
2952 allocate(sigmap1_scbase(ntyp_molec(1),ntyp_molec(2)))
2953 allocate(sigmap2_scbase(ntyp_molec(1),ntyp_molec(2)))
2954 allocate(chis_scbase(ntyp_molec(1),ntyp_molec(2),2))
2955 allocate(dhead_scbasei(ntyp_molec(1),ntyp_molec(2)))
2956 allocate(dhead_scbasej(ntyp_molec(1),ntyp_molec(2)))
2957 allocate(rborn_scbasei(ntyp_molec(1),ntyp_molec(2)))
2958 allocate(rborn_scbasej(ntyp_molec(1),ntyp_molec(2)))
2959 allocate(wdipdip_scbase(3,ntyp_molec(1),ntyp_molec(2)))
2960 allocate(wqdip_scbase(2,ntyp_molec(1),ntyp_molec(2)))
2961 allocate(alphapol_scbase(ntyp_molec(1),ntyp_molec(2)))
2962 allocate(epsintab_scbase(ntyp_molec(1),ntyp_molec(2)))
2965 do i=1,ntyp_molec(1)
2966 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
2967 write (*,*) "Im in ", i, " ", j
2968 read(isidep_scbase,*) &
2969 eps_scbase(i,j),sigma_scbase(i,j),chi_scbase(i,j,1),&
2970 chi_scbase(i,j,2),chipp_scbase(i,j,1),chipp_scbase(i,j,2)
2971 write(*,*) "eps",eps_scbase(i,j)
2972 read(isidep_scbase,*) &
2973 (alphasur_scbase(k,i,j),k=1,4),sigmap1_scbase(i,j),sigmap2_scbase(i,j), &
2974 chis_scbase(i,j,1),chis_scbase(i,j,2)
2975 read(isidep_scbase,*) &
2976 dhead_scbasei(i,j), &
2977 dhead_scbasej(i,j), &
2978 rborn_scbasei(i,j),rborn_scbasej(i,j)
2979 read(isidep_scbase,*) &
2980 (wdipdip_scbase(k,i,j),k=1,3), &
2981 (wqdip_scbase(k,i,j),k=1,2)
2982 read(isidep_scbase,*) &
2983 alphapol_scbase(i,j), &
2984 epsintab_scbase(i,j)
2987 allocate(aa_scbase(ntyp_molec(1),ntyp_molec(2)))
2988 allocate(bb_scbase(ntyp_molec(1),ntyp_molec(2)))
2990 do i=1,ntyp_molec(1)
2991 do j=1,ntyp_molec(2)-1
2992 epsij=eps_scbase(i,j)
2993 rrij=sigma_scbase(i,j)
2998 sigeps=dsign(1.0D0,epsij)
3000 aa_scbase(i,j)=epsij*rrij*rrij
3001 bb_scbase(i,j)=-sigeps*epsij*rrij
3004 !-----------------READING PEP BASE POTENTIALS-------------------
3005 allocate(eps_pepbase(ntyp_molec(2)))
3006 allocate(sigma_pepbase(ntyp_molec(2)))
3007 allocate(chi_pepbase(ntyp_molec(2),2))
3008 allocate(chipp_pepbase(ntyp_molec(2),2))
3009 allocate(alphasur_pepbase(4,ntyp_molec(2)))
3010 allocate(sigmap1_pepbase(ntyp_molec(2)))
3011 allocate(sigmap2_pepbase(ntyp_molec(2)))
3012 allocate(chis_pepbase(ntyp_molec(2),2))
3013 allocate(wdipdip_pepbase(3,ntyp_molec(2)))
3016 do j=1,ntyp_molec(2)-1 ! without U then we will take T for U
3017 write (*,*) "Im in ", i, " ", j
3018 read(isidep_pepbase,*) &
3019 eps_pepbase(j),sigma_pepbase(j),chi_pepbase(j,1),&
3020 chi_pepbase(j,2),chipp_pepbase(j,1),chipp_pepbase(j,2)
3021 write(*,*) "eps",eps_pepbase(j)
3022 read(isidep_pepbase,*) &
3023 (alphasur_pepbase(k,j),k=1,4),sigmap1_pepbase(j),sigmap2_pepbase(j), &
3024 chis_pepbase(j,1),chis_pepbase(j,2)
3025 read(isidep_pepbase,*) &
3026 (wdipdip_pepbase(k,j),k=1,3)
3028 allocate(aa_pepbase(ntyp_molec(2)))
3029 allocate(bb_pepbase(ntyp_molec(2)))
3031 do j=1,ntyp_molec(2)-1
3032 epsij=eps_pepbase(j)
3033 rrij=sigma_pepbase(j)
3038 sigeps=dsign(1.0D0,epsij)
3040 aa_pepbase(j)=epsij*rrij*rrij
3041 bb_pepbase(j)=-sigeps*epsij*rrij
3043 !--------------READING SC PHOSPHATE-------------------------------------
3044 allocate(eps_scpho(ntyp_molec(1)))
3045 allocate(sigma_scpho(ntyp_molec(1)))
3046 allocate(chi_scpho(ntyp_molec(1),2))
3047 allocate(chipp_scpho(ntyp_molec(1),2))
3048 allocate(alphasur_scpho(4,ntyp_molec(1)))
3049 allocate(sigmap1_scpho(ntyp_molec(1)))
3050 allocate(sigmap2_scpho(ntyp_molec(1)))
3051 allocate(chis_scpho(ntyp_molec(1),2))
3052 allocate(wqq_scpho(ntyp_molec(1)))
3053 allocate(wqdip_scpho(2,ntyp_molec(1)))
3054 allocate(alphapol_scpho(ntyp_molec(1)))
3055 allocate(epsintab_scpho(ntyp_molec(1)))
3056 allocate(dhead_scphoi(ntyp_molec(1)))
3057 allocate(rborn_scphoi(ntyp_molec(1)))
3058 allocate(rborn_scphoj(ntyp_molec(1)))
3059 allocate(alphi_scpho(ntyp_molec(1)))
3063 do j=1,ntyp_molec(1) ! without U then we will take T for U
3064 write (*,*) "Im in scpho ", i, " ", j
3065 read(isidep_scpho,*) &
3066 eps_scpho(j),sigma_scpho(j),chi_scpho(j,1),&
3067 chi_scpho(j,2),chipp_scpho(j,1),chipp_scpho(j,2)
3068 write(*,*) "eps",eps_scpho(j)
3069 read(isidep_scpho,*) &
3070 (alphasur_scpho(k,j),k=1,4),sigmap1_scpho(j),sigmap2_scpho(j), &
3071 chis_scpho(j,1),chis_scpho(j,2)
3072 read(isidep_scpho,*) &
3073 (wqdip_scpho(k,j),k=1,2),wqq_scpho(j),dhead_scphoi(j)
3074 read(isidep_scpho,*) &
3075 epsintab_scpho(j),alphapol_scpho(j),rborn_scphoi(j),rborn_scphoj(j), &
3079 allocate(aa_scpho(ntyp_molec(1)))
3080 allocate(bb_scpho(ntyp_molec(1)))
3082 do j=1,ntyp_molec(1)
3089 sigeps=dsign(1.0D0,epsij)
3091 aa_scpho(j)=epsij*rrij*rrij
3092 bb_scpho(j)=-sigeps*epsij*rrij
3096 read(isidep_peppho,*) &
3097 eps_peppho,sigma_peppho
3098 read(isidep_peppho,*) &
3099 (alphasur_peppho(k),k=1,4),sigmap1_peppho,sigmap2_peppho
3100 read(isidep_peppho,*) &
3101 (wqdip_peppho(k),k=1,2)
3109 sigeps=dsign(1.0D0,epsij)
3111 aa_peppho=epsij*rrij*rrij
3112 bb_peppho=-sigeps*epsij*rrij
3115 allocate(aad(ntyp,2),bad(ntyp,2)) !(ntyp,2)
3120 ! Define the SC-p interaction constants (hard-coded; old style)
3123 ! "Soft" SC-p repulsion (causes helices to be too flat, but facilitates
3125 ! aad(i,1)=0.3D0*4.0D0**12
3126 ! Following line for constants currently implemented
3127 ! "Hard" SC-p repulsion (gives correct turn spacing in helices)
3128 aad(i,1)=1.5D0*4.0D0**12
3129 ! aad(i,1)=0.17D0*5.6D0**12
3131 ! "Soft" SC-p repulsion
3133 ! Following line for constants currently implemented
3134 ! aad(i,1)=0.3D0*4.0D0**6
3135 ! "Hard" SC-p repulsion
3136 bad(i,1)=3.0D0*4.0D0**6
3137 ! bad(i,1)=-2.0D0*0.17D0*5.6D0**6
3146 ! 8/9/01 Read the SC-p interaction constants from file
3149 read (iscpp,*,end=118,err=118) (eps_scp(i,j),rscp(i,j),j=1,2)
3152 aad(i,1)=dabs(eps_scp(i,1))*rscp(i,1)**12
3153 aad(i,2)=dabs(eps_scp(i,2))*rscp(i,2)**12
3154 bad(i,1)=-2*eps_scp(i,1)*rscp(i,1)**6
3155 bad(i,2)=-2*eps_scp(i,2)*rscp(i,2)**6
3159 write (iout,*) "Parameters of SC-p interactions:"
3161 write (iout,'(4f8.3,4e12.4)') eps_scp(i,1),rscp(i,1),&
3162 eps_scp(i,2),rscp(i,2),aad(i,1),bad(i,1),aad(i,2),bad(i,2)
3167 allocate(aad_nucl(ntyp_molec(2)),bad_nucl(ntyp_molec(2))) !(ntyp,2)
3169 do i=1,ntyp_molec(2)
3170 read (iscpp_nucl,*,end=118,err=118) eps_scp_nucl(i),rscp_nucl(i)
3172 do i=1,ntyp_molec(2)
3173 aad_nucl(i)=dabs(eps_scp_nucl(i))*rscp_nucl(i)**12
3174 bad_nucl(i)=-2*eps_scp_nucl(i)*rscp_nucl(i)**6
3176 r0pp=1.12246204830937298142*5.16158
3182 ! Define the constants of the disulfide bridge
3186 ! Old arbitrary potential - commented out.
3191 ! Constants of the disulfide-bond potential determined based on the RHF/6-31G**
3192 ! energy surface of diethyl disulfide.
3193 ! A. Liwo and U. Kozlowska, 11/24/03
3210 write (iout,'(/a)') "Disulfide bridge parameters:"
3211 write (iout,'(a,f10.2)') 'S-S bridge energy: ',ebr
3212 write (iout,'(2(a,f10.2))') 'd0cm:',d0cm,' akcm:',akcm
3213 write (iout,'(2(a,f10.2))') 'akth:',akth,' akct:',akct
3214 write (iout,'(3(a,f10.2))') 'v1ss:',v1ss,' v2ss:',v2ss,&
3217 if (shield_mode.gt.0) then
3218 pi=4.0D0*datan(1.0D0)
3219 !C VSolvSphere the volume of solving sphere
3221 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
3222 !C there will be no distinction between proline peptide group and normal peptide
3223 !C group in case of shielding parameters
3224 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
3225 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
3226 write (iout,*) VSolvSphere,VSolvSphere_div
3227 !C long axis of side chain
3229 long_r_sidechain(i)=vbldsc0(1,i)
3230 ! if (scelemode.eq.0) then
3231 short_r_sidechain(i)=sigma(i,i)/sqrt(2.0)
3232 if (short_r_sidechain(i).eq.0.0) short_r_sidechain(i)=0.2
3234 ! short_r_sidechain(i)=sigma(i,i)
3236 write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
3243 111 write (iout,*) "Error reading bending energy parameters."
3245 112 write (iout,*) "Error reading rotamer energy parameters."
3247 113 write (iout,*) "Error reading torsional energy parameters."
3249 114 write (iout,*) "Error reading double torsional energy parameters."
3251 115 write (iout,*) &
3252 "Error reading cumulant (multibody energy) parameters."
3254 116 write (iout,*) "Error reading electrostatic energy parameters."
3256 117 write (iout,*) "Error reading side chain interaction parameters."
3258 118 write (iout,*) "Error reading SCp interaction parameters."
3260 119 write (iout,*) "Error reading SCCOR parameters"
3262 121 write (iout,*) "Error in Czybyshev parameters"
3265 call MPI_Finalize(Ierror)
3269 end subroutine parmread
3271 !-----------------------------------------------------------------------------
3273 !-----------------------------------------------------------------------------
3274 subroutine printmat(ldim,m,n,iout,key,a)
3277 character(len=3),dimension(n) :: key
3278 real(kind=8),dimension(ldim,n) :: a
3280 integer :: i,j,k,m,iout,nlim
3284 write (iout,1000) (key(k),k=i,nlim)
3286 1000 format (/5x,8(6x,a3))
3287 1020 format (/80(1h-)/)
3289 write (iout,1010) key(j),(a(j,k),k=i,nlim)
3292 1010 format (a3,2x,8(f9.4))
3294 end subroutine printmat
3295 !-----------------------------------------------------------------------------
3297 !-----------------------------------------------------------------------------
3299 ! Read the PDB file and convert the peptide geometry into virtual-chain
3302 use energy_data, only: itype,istype
3306 ! use control, only: rescode,sugarcode
3307 ! implicit real*8 (a-h,o-z)
3308 ! include 'DIMENSIONS'
3309 ! include 'COMMON.LOCAL'
3310 ! include 'COMMON.VAR'
3311 ! include 'COMMON.CHAIN'
3312 ! include 'COMMON.INTERACT'
3313 ! include 'COMMON.IOUNITS'
3314 ! include 'COMMON.GEO'
3315 ! include 'COMMON.NAMES'
3316 ! include 'COMMON.CONTROL'
3317 ! include 'COMMON.DISTFIT'
3318 ! include 'COMMON.SETUP'
3319 integer :: i,j,ibeg,ishift1,ires,iii,ires_old,ishift,k!,ity!,&
3321 logical :: lprn=.true.,fail
3322 real(kind=8),dimension(3) :: e1,e2,e3
3323 real(kind=8) :: dcj,efree_temp
3324 character(len=3) :: seq,res,res2
3325 character(len=5) :: atom
3326 character(len=80) :: card
3327 real(kind=8),dimension(3,20) :: sccor
3328 integer :: kkk,lll,icha,kupa,molecule,counter,seqalingbegin !rescode,
3329 integer :: isugar,molecprev,firstion
3330 character*1 :: sugar
3332 real(kind=8),dimension(3) :: ccc
3334 integer,dimension(2,maxres/3) :: hfrag_alloc
3335 integer,dimension(4,maxres/3) :: bfrag_alloc
3336 real(kind=8),dimension(3,maxres2+2,maxperm) :: cref_alloc !(3,maxres2+2,maxperm)
3337 real(kind=8),dimension(:,:), allocatable :: c_temporary
3338 integer,dimension(:,:) , allocatable :: itype_temporary
3339 integer,dimension(:),allocatable :: istype_temp
3346 ! write (2,*) "UNRES_PDB",unres_pdb
3366 !-----------------------------
3367 allocate(hfrag(2,maxres/3)) !(2,maxres/3)
3368 allocate(bfrag(4,maxres/3)) !(4,maxres/3)
3369 if(.not. allocated(istype)) allocate(istype(maxres))
3371 read (ipdbin,'(a80)',end=10) card
3372 write (iout,'(a)') card
3373 if (card(:5).eq.'HELIX') then
3376 read(card(22:25),*) hfrag(1,nhfrag)
3377 read(card(34:37),*) hfrag(2,nhfrag)
3379 if (card(:5).eq.'SHEET') then
3382 read(card(24:26),*) bfrag(1,nbfrag)
3383 read(card(35:37),*) bfrag(2,nbfrag)
3384 !rc----------------------------------------
3385 !rc to be corrected !!!
3386 bfrag(3,nbfrag)=bfrag(1,nbfrag)
3387 bfrag(4,nbfrag)=bfrag(2,nbfrag)
3388 !rc----------------------------------------
3390 if (card(:3).eq.'END') then
3392 else if (card(:3).eq.'TER') then
3397 itype(ires_old,molecule)=ntyp1_molec(molecule)
3398 itype(ires_old-1,molecule)=ntyp1_molec(molecule)
3399 nres_molec(molecule)=nres_molec(molecule)+2
3401 ! write (iout,*) "Chain ended",ires,ishift,ires_old
3404 dc(j,ires)=sccor(j,iii)
3407 call sccenter(ires,iii,sccor)
3413 if (index(card,"FREE ENERGY").gt.0) read(card(35:),*) efree_temp
3414 ! Fish out the ATOM cards.
3415 ! write(iout,*) 'card',card(1:20)
3416 ! print *,"ATU ",card(1:6), CARD(3:6)
3418 if (index(card(1:4),'ATOM').gt.0) then
3419 read (card(12:16),*) atom
3420 ! write (iout,*) "! ",atom," !",ires
3421 ! if (atom.eq.'CA' .or. atom.eq.'CH3') then
3422 read (card(23:26),*) ires
3423 read (card(18:20),'(a3)') res
3424 ! write (iout,*) "ires",ires,ires-ishift+ishift1,
3425 ! & " ires_old",ires_old
3426 ! write (iout,*) "ishift",ishift," ishift1",ishift1
3427 ! write (iout,*) "IRES",ires-ishift+ishift1,ires_old
3428 if (ires-ishift+ishift1.ne.ires_old) then
3429 ! Calculate the CM of the preceding residue.
3430 ! if (ibeg.eq.0) call sccenter(ires,iii,sccor)
3432 ! write (iout,*) "Calculating sidechain center iii",iii
3435 dc(j,ires+ishift1-ishift-1)=sccor(j,iii)
3438 call sccenter(ires_old,iii,sccor)
3442 ! Start new residue.
3443 if (res.eq.'Cl-' .or. res.eq.'Na+') then
3446 else if (ibeg.eq.1) then
3447 write (iout,*) "BEG ires",ires
3449 if (res.ne.'GLY' .and. res.ne. 'ACE') then
3452 nres_molec(molecule)=nres_molec(molecule)+1
3454 ires=ires-ishift+ishift1
3456 ! write (iout,*) "ishift",ishift," ires",ires,&
3457 ! " ires_old",ires_old
3459 else if (ibeg.eq.2) then
3461 ishift=-ires_old+ires-1 !!!!!
3462 ishift1=ishift1-1 !!!!!
3463 ! write (iout,*) "New chain started",ires,ishift,ishift1,"!"
3464 ires=ires-ishift+ishift1
3465 ! print *,ires,ishift,ishift1
3469 ishift=ishift-(ires-ishift+ishift1-ires_old-1)
3470 ires=ires-ishift+ishift1
3473 ! print *,'atom',ires,atom
3474 if (res.eq.'ACE' .or. res.eq.'NHE') then
3477 if (atom.eq.'CA '.or.atom.eq.'N ') then
3479 itype(ires,molecule)=rescode(ires,res,0,molecule)
3481 ! nres_molec(molecule)=nres_molec(molecule)+1
3485 itype(ires,molecule)=rescode(ires,res2,0,molecule)
3486 ! nres_molec(molecule)=nres_molec(molecule)+1
3487 read (card(19:19),'(a1)') sugar
3488 isugar=sugarcode(sugar,ires)
3489 ! if (ibeg.eq.1) then
3493 ! print *,"ires=",ires,istype(ires)
3499 ires=ires-ishift+ishift1
3501 ! write (iout,*) "ires_old",ires_old," ires",ires
3502 if (card(27:27).eq."A" .or. card(27:27).eq."B") then
3505 ! write (2,*) "ires",ires," res ",res!," ity"!,ity
3506 if (atom.eq.'CA' .or. atom.eq.'CH3' .or. &
3507 res.eq.'NHE'.and.atom(:2).eq.'HN') then
3508 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3509 ! print *,ires,ishift,ishift1
3510 ! write (iout,*) "backbone ",atom
3512 write (iout,'(2i3,2x,a,3f8.3)') &
3513 ires,itype(ires,1),res,(c(j,ires),j=1,3)
3516 nres_molec(molecule)=nres_molec(molecule)+1
3518 sccor(j,iii)=c(j,ires)
3520 else if (.not.unres_pdb .and. (atom.eq."C1'" .or. &
3521 atom.eq."C2'" .or. atom.eq."C3'" &
3522 .or. atom.eq."C4'" .or. atom.eq."O4'")) then
3523 read(card(31:54),'(3f8.3)') (ccc(j),j=1,3)
3524 !c write (2,'(i5,3f10.5)') ires,(ccc(j),j=1,3)
3525 ! print *,ires,ishift,ishift1
3529 ! sccor(j,iii)=c(j,ires)
3532 c(j,ires)=c(j,ires)+ccc(j)/5.0
3534 print *,counter,molecule
3535 if (counter.eq.5) then
3537 nres_molec(molecule)=nres_molec(molecule)+1
3540 ! sccor(j,iii)=c(j,ires)
3544 ! print *, "ATOM",atom(1:3)
3545 else if (atom.eq."C5'") then
3546 read (card(19:19),'(a1)') sugar
3547 isugar=sugarcode(sugar,ires)
3552 ! print *,ires,istype(ires)
3556 ! print *,"nres_molec(molecule)",nres_molec(molecule),ires
3557 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3558 nres_molec(molecule)=nres_molec(molecule)+1
3559 print *,"nres_molec(molecule)",nres_molec(molecule),ires
3563 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3565 else if ((atom.eq."C1'").and.unres_pdb) then
3567 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3568 ! write (*,*) card(23:27),ires,itype(ires,1)
3569 else if (atom.ne.'O'.and.atom(1:1).ne.'H' .and. &
3570 atom.ne.'N' .and. atom.ne.'C' .and. &
3571 atom(:2).ne.'1H' .and. atom(:2).ne.'2H' .and. &
3572 atom.ne.'OXT' .and. atom(:2).ne.'3H' &
3573 .and. atom.ne.'P '.and. &
3574 atom(1:1).ne.'H' .and. &
3575 atom.ne.'OP1' .and. atom.ne.'OP2 '.and. atom.ne.'OP3'&
3577 ! write (iout,*) "sidechain ",atom
3578 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3579 if ((molecule.ne.2).or.(atom(3:3).ne."'")) then
3580 ! write (iout,*) "sidechain ",atom,molecule,ires,atom(3:3)
3582 read(card(31:54),'(3f8.3)') (sccor(j,iii),j=1,3)
3585 ! print *,"IONS",ions,card(1:6)
3586 else if ((ions).and.(card(1:6).eq.'HETATM')) then
3587 if (firstion.eq.0) then
3591 dc(j,ires)=sccor(j,iii)
3594 call sccenter(ires,iii,sccor)
3597 read (card(12:16),*) atom
3598 ! print *,"HETATOM", atom
3599 read (card(18:20),'(a3)') res
3600 if ((atom(1:2).eq.'NA').or.(atom(1:2).eq.'CL').or.&
3601 (atom(1:2).eq.'CA').or.(atom(1:2).eq.'MG') &
3602 .or.(atom(1:2).eq.'K ')) &
3605 if (molecule.ne.5) molecprev=molecule
3607 nres_molec(molecule)=nres_molec(molecule)+1
3608 print *,"HERE",nres_molec(molecule)
3610 itype(ires,molecule)=rescode(ires,res,0,molecule)
3611 read(card(31:54),'(3f8.3)') (c(j,ires),j=1,3)
3615 10 write (iout,'(a,i5)') ' Number of residues found: ',ires
3616 if (ires.eq.0) return
3617 ! Calculate dummy residue coordinates inside the "chain" of a multichain
3620 if (((ires_old.ne.ires).and.(molecule.ne.5)) &
3622 nres_molec(molecule)=nres_molec(molecule)-2
3623 print *,'I have',nres, nres_molec(:)
3625 do k=1,4 ! ions are without dummy
3626 if (nres_molec(k).eq.0) cycle
3628 ! write (iout,*) i,itype(i,1)
3629 ! if (itype(i,1).eq.ntyp1) then
3630 ! write (iout,*) "dummy",i,itype(i,1)
3632 ! c(j,i)=((c(j,i-1)+c(j,i+1))/2+2*c(j,i-1)-c(j,i-2))/2
3633 ! c(j,i)=(c(j,i-1)+c(j,i+1))/2
3637 if (itype(i,k).eq.ntyp1_molec(k)) then
3638 if (itype(i+1,k).eq.ntyp1_molec(k)) then
3639 if (itype(i+2,k).eq.0) then
3640 ! print *,"masz sieczke"
3642 if (itype(i+2,j).ne.0) then
3644 itype(i+1,j)=ntyp1_molec(j)
3645 nres_molec(k)=nres_molec(k)-1
3646 nres_molec(j)=nres_molec(j)+1
3652 ! 16/01/2014 by Adasko: Adding to dummy atoms in the chain
3653 ! first is connected prevous chain (itype(i+1,1).eq.ntyp1)=true
3654 ! second dummy atom is conected to next chain itype(i+1,1).eq.ntyp1=false
3655 ! if (unres_pdb) then
3656 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3657 ! print *,i,'tu dochodze'
3658 ! call refsys(i-3,i-2,i-1,e1,e2,e3,fail)
3666 ! c(j,i)=c(j,i-1)-1.9d0*e2(j)
3670 dcj=(c(j,i-2)-c(j,i-3))/2.0
3671 if (dcj.eq.0) dcj=1.23591524223
3676 else !itype(i+1,1).eq.ntyp1
3677 ! if (unres_pdb) then
3678 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3679 ! call refsys(i+1,i+2,i+3,e1,e2,e3,fail)
3686 ! c(j,i)=c(j,i+1)-1.9d0*e2(j)
3687 c(j,i)=c(j,i-1)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
3691 dcj=(c(j,i+3)-c(j,i+2))/2.0
3692 if (dcj.eq.0) dcj=1.23591524223
3697 endif !itype(i+1,1).eq.ntyp1
3698 endif !itype.eq.ntyp1
3702 ! Calculate the CM of the last side chain.
3706 dc(j,ires)=sccor(j,iii)
3709 call sccenter(ires,iii,sccor)
3715 ! print *,"molecule",molecule
3716 if ((itype(nres,1).ne.10)) then
3718 if (molecule.eq.5) molecule=molecprev
3719 itype(nres,molecule)=ntyp1_molec(molecule)
3720 nres_molec(molecule)=nres_molec(molecule)+1
3721 ! if (unres_pdb) then
3722 ! 2/15/2013 by Adam: corrected insertion of the last dummy residue
3723 ! call refsys(nres-3,nres-2,nres-1,e1,e2,e3,fail)
3730 ! c(j,nres)=c(j,nres-1)-1.9d0*e2(j)
3734 dcj=(c(j,nres-2)-c(j,nres-3))/2.0
3735 c(j,nres)=c(j,nres-1)+dcj
3736 c(j,2*nres)=c(j,nres)
3740 ! print *,'I have',nres, nres_molec(:)
3742 !el kontrola nres w pliku inputowym WHAM-a w porownaniu z wartoscia wczytana z pliku pdb
3744 if (nres.ne.nres0) then
3745 write (iout,*) "Error: wrong parameter value: NRES=",nres,&
3747 stop "Error nres value in WHAM input"
3750 !---------------------------------
3751 !el reallocate tables
3754 ! hfrag_alloc(j,i)=hfrag(j,i)
3757 ! bfrag_alloc(j,i)=bfrag(j,i)
3763 ! allocate(hfrag(2,nres/3)) !(2,maxres/3)
3764 !el allocate(hfrag(2,nhfrag)) !(2,maxres/3)
3765 !el allocate(bfrag(4,nbfrag)) !(4,maxres/3)
3766 ! allocate(bfrag(4,nres/3)) !(4,maxres/3)
3770 ! hfrag(j,i)=hfrag_alloc(j,i)
3775 ! bfrag(j,i)=bfrag_alloc(j,i)
3778 !el end reallocate tables
3779 !---------------------------------
3787 c(j,2*nres)=c(j,nres)
3790 if (itype(1,1).eq.ntyp1) then
3794 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3795 call refsys(2,3,4,e1,e2,e3,fail)
3802 ! c(j,1)=c(j,2)-1.9d0*e2(j)
3803 c(j,1)=c(j,2)+1.9d0*(-e1(j)+e2(j))/sqrt(2.0d0)
3807 dcj=(c(j,4)-c(j,3))/2.0
3813 ! First lets assign correct dummy to correct type of chain
3815 if (itype(1,1).eq.ntyp1) then
3816 if (itype(2,1).eq.0) then
3818 if (itype(2,j).ne.0) then
3820 itype(1,j)=ntyp1_molec(j)
3821 nres_molec(1)=nres_molec(1)-1
3822 nres_molec(j)=nres_molec(j)+1
3829 print *,'I have',nres, nres_molec(:)
3831 ! Copy the coordinates to reference coordinates
3837 ! Calculate internal coordinates.
3839 write (iout,'(/a)') &
3840 "Cartesian coordinates of the reference structure"
3841 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
3842 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3844 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
3845 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3846 (c(j,ires+nres),j=1,3)
3849 ! znamy już nres wiec mozna alokowac tablice
3850 ! Calculate internal coordinates.
3851 if(me.eq.king.or..not.out1file)then
3852 write (iout,'(a)') &
3853 "Backbone and SC coordinates as read from the PDB"
3855 write (iout,'(i5,i3,2x,a,3f8.3,5x,3f8.3)') &
3856 ires,itype(ires,1),restyp(itype(ires,1),1),(c(j,ires),j=1,3),&
3857 (c(j,nres+ires),j=1,3)
3860 ! NOW LETS ROCK! SORTING
3861 allocate(c_temporary(3,2*nres))
3862 allocate(itype_temporary(nres,5))
3863 if (.not.allocated(molnum)) allocate(molnum(nres+1))
3864 if (.not.allocated(istype)) write(iout,*) &
3865 "SOMETHING WRONG WITH ISTYTPE"
3866 allocate(istype_temp(nres))
3867 itype_temporary(:,:)=0
3871 if (itype(i,k).ne.0) then
3873 c_temporary(j,seqalingbegin)=c(j,i)
3874 c_temporary(j,seqalingbegin+nres)=c(j,i+nres)
3877 itype_temporary(seqalingbegin,k)=itype(i,k)
3878 print *,i,k,itype(i,k),itype_temporary(seqalingbegin,k),seqalingbegin
3879 istype_temp(seqalingbegin)=istype(i)
3880 molnum(seqalingbegin)=k
3881 seqalingbegin=seqalingbegin+1
3887 c(j,i)=c_temporary(j,i)
3892 itype(i,k)=itype_temporary(i,k)
3893 istype(i)=istype_temp(i)
3896 ! if (itype(1,1).eq.ntyp1) then
3899 ! if (unres_pdb) then
3900 ! 2/15/2013 by Adam: corrected insertion of the first dummy residue
3901 ! call refsys(2,3,4,e1,e2,e3,fail)
3908 ! c(j,1)=c(j,2)-1.9d0*e2(j)
3912 ! dcj=(c(j,4)-c(j,3))/2.0
3914 ! c(j,nres+1)=c(j,1)
3920 write (iout,'(/a)') &
3921 "Cartesian coordinates of the reference structure after sorting"
3922 write (iout,'(a,16x,3(3x,a5),5x,3(3x,a5))') &
3923 "Residue","X(CA)","Y(CA)","Z(CA)","X(SC)","Y(SC)","Z(SC)"
3925 write (iout,'(5(a3,1x),i5,3f8.3,5x,3f8.3)') &
3926 (restyp(itype(ires,j),j),j=1,5),ires,(c(j,ires),j=1,3),&
3927 (c(j,ires+nres),j=1,3)
3931 ! print *,seqalingbegin,nres
3932 if(.not.allocated(vbld)) then
3933 allocate(vbld(2*nres))
3938 if(.not.allocated(vbld_inv)) then
3939 allocate(vbld_inv(2*nres))
3945 if(.not.allocated(theta)) then
3946 allocate(theta(nres+2))
3950 if(.not.allocated(phi)) allocate(phi(nres+2))
3951 if(.not.allocated(alph)) allocate(alph(nres+2))
3952 if(.not.allocated(omeg)) allocate(omeg(nres+2))
3953 if(.not.allocated(thetaref)) allocate(thetaref(nres+2))
3954 if(.not.allocated(phiref)) allocate(phiref(nres+2))
3955 if(.not.allocated(costtab)) allocate(costtab(nres))
3956 if(.not.allocated(sinttab)) allocate(sinttab(nres))
3957 if(.not.allocated(cost2tab)) allocate(cost2tab(nres))
3958 if(.not.allocated(sint2tab)) allocate(sint2tab(nres))
3959 if(.not.allocated(xxref)) allocate(xxref(nres))
3960 if(.not.allocated(yyref)) allocate(yyref(nres))
3961 if(.not.allocated(zzref)) allocate(zzref(nres)) !(maxres)
3962 if(.not.allocated(dc_norm)) then
3963 ! if(.not.allocated(dc_norm)) allocate(dc_norm(3,0:2*nres+2))
3964 allocate(dc_norm(3,0:2*nres+2))
3968 call int_from_cart(.true.,.false.)
3969 call sc_loc_geom(.false.)
3971 thetaref(i)=theta(i)
3981 dc(j,i)=c(j,i+1)-c(j,i)
3982 dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
3987 dc(j,i+nres)=c(j,i+nres)-c(j,i)
3988 dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
3990 ! write (iout,*) i,(dc(j,i+nres),j=1,3),(dc_norm(j,i+nres),j=1,3),&
3994 ! Copy the coordinates to reference coordinates
3995 ! Splits to single chain if occurs
3999 ! cref(j,i,cou)=c(j,i)
4003 if(.not.allocated(cref)) allocate(cref(3,2*nres+2,maxperm)) !(3,maxres2+2,maxperm)
4004 if(.not.allocated(chain_rep)) allocate(chain_rep(3,2*nres+2,maxsym)) !(3,maxres2+2,maxsym)
4005 if(.not.allocated(tabperm)) allocate(tabperm(maxperm,maxsym)) !(maxperm,maxsym)
4006 !-----------------------------
4010 write (iout,*) "symetr", symetr
4013 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4015 if ((itype(i-1,1).eq.ntyp1).and.(i.gt.2)) then
4018 ! write (iout,*) "spraw lancuchy",(c(j,i),j=1,3)
4023 cref(j,i,cou)=c(j,i)
4024 cref(j,i+nres,cou)=c(j,i+nres)
4026 chain_rep(j,lll,kkk)=c(j,i)
4027 chain_rep(j,lll+nres,kkk)=c(j,i+nres)
4031 write (iout,*) chain_length
4032 if (chain_length.eq.0) chain_length=nres
4034 chain_rep(j,chain_length,symetr)=chain_rep(j,chain_length,1)
4035 chain_rep(j,chain_length+nres,symetr) &
4036 =chain_rep(j,chain_length+nres,1)
4039 ! write (iout,*) "spraw lancuchy",chain_length,symetr
4041 ! do kkk=1,chain_length
4042 ! write (iout,*) itype(kkk,1),(chain_rep(j,kkk,i), j=1,3)
4046 ! makes copy of chains
4047 write (iout,*) "symetr", symetr
4052 if (symetr.gt.1) then
4059 write(iout,*) (tabperm(i,kkk),kkk=1,4)
4065 write (iout,*) i,icha
4066 do lll=1,chain_length
4068 if (cou.le.nres) then
4070 kupa=mod(lll,chain_length)
4071 iprzes=(kkk-1)*chain_length+lll
4072 if (kupa.eq.0) kupa=chain_length
4073 write (iout,*) "kupa", kupa
4074 cref(j,iprzes,i)=chain_rep(j,kupa,icha)
4075 cref(j,iprzes+nres,i)=chain_rep(j,kupa+nres,icha)
4082 !-koniec robienia kopii
4085 write (iout,*) "nowa struktura", nperm
4087 write (iout,110) restyp(itype(i,1),1),i,cref(1,i,kkk),&
4089 cref(3,i,kkk),cref(1,nres+i,kkk),&
4090 cref(2,nres+i,kkk),cref(3,nres+i,kkk)
4092 100 format (//' alpha-carbon coordinates ',&
4093 ' centroid coordinates'/ &
4094 ' ', 6X,'X',11X,'Y',11X,'Z', &
4095 10X,'X',11X,'Y',11X,'Z')
4096 110 format (a,'(',i5,')',6f12.5)
4102 bfrag(i,j)=bfrag(i,j)-ishift
4108 hfrag(i,j)=hfrag(i,j)-ishift
4114 end subroutine readpdb
4115 #if .not. defined(WHAM_RUN) && .not. defined(CLUSTER)
4116 !-----------------------------------------------------------------------------
4118 !-----------------------------------------------------------------------------
4119 subroutine read_control
4133 use random, only: random_init
4134 ! implicit real*8 (a-h,o-z)
4135 ! include 'DIMENSIONS'
4137 use prng, only:prng_restart
4139 logical :: OKRandom!, prng_restart
4142 ! include 'COMMON.IOUNITS'
4143 ! include 'COMMON.TIME1'
4144 ! include 'COMMON.THREAD'
4145 ! include 'COMMON.SBRIDGE'
4146 ! include 'COMMON.CONTROL'
4147 ! include 'COMMON.MCM'
4148 ! include 'COMMON.MAP'
4149 ! include 'COMMON.HEADER'
4150 ! include 'COMMON.CSA'
4151 ! include 'COMMON.CHAIN'
4152 ! include 'COMMON.MUCA'
4153 ! include 'COMMON.MD'
4154 ! include 'COMMON.FFIELD'
4155 ! include 'COMMON.INTERACT'
4156 ! include 'COMMON.SETUP'
4157 !el integer :: KDIAG,ICORFL,IXDR
4158 !el COMMON /MACHSW/ KDIAG,ICORFL,IXDR
4159 character(len=8),dimension(0:3) :: diagmeth = reshape((/'Library ',&
4160 'EVVRSP ','Givens ','Jacobi '/),shape(diagmeth))
4161 ! character(len=80) :: ucase
4162 character(len=640) :: controlcard
4164 real(kind=8) :: seed,rmsdbc,rmsdbc1max,rmsdbcm,drms,timem!,&
4170 read (INP,'(a)') titel
4171 call card_concat(controlcard,.true.)
4172 ! out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
4173 ! print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
4174 call reada(controlcard,'SEED',seed,0.0D0)
4175 call random_init(seed)
4176 ! Set up the time limit (caution! The time must be input in minutes!)
4177 read_cart=index(controlcard,'READ_CART').gt.0
4178 call readi(controlcard,'CONSTR_DIST',constr_dist,0)
4179 call readi(controlcard,'SYM',symetr,1)
4180 call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
4181 unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
4182 call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
4183 call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
4184 call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
4185 call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
4186 call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
4187 call reada(controlcard,'DRMS',drms,0.1D0)
4188 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4189 write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
4190 write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
4191 write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
4192 write (iout,'(a,f10.1)')'DRMS = ',drms
4193 write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
4194 write (iout,'(a,f10.1)') 'Time limit (min):',timlim
4196 call readi(controlcard,'NZ_START',nz_start,0)
4197 call readi(controlcard,'NZ_END',nz_end,0)
4198 call readi(controlcard,'IZ_SC',iz_sc,0)
4199 timlim=60.0D0*timlim
4200 safety = 60.0d0*safety
4203 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4204 !C SHIELD keyword sets if the shielding effect of side-chains is used
4205 !C 0 denots no shielding is used all peptide are equally despite the
4206 !C solvent accesible area
4207 !C 1 the newly introduced function
4208 !C 2 reseved for further possible developement
4209 call readi(controlcard,'SHIELD',shield_mode,0)
4210 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4211 write(iout,*) "shield_mode",shield_mode
4212 call readi(controlcard,'TORMODE',tor_mode,0)
4213 !C if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
4214 write(iout,*) "torsional and valence angle mode",tor_mode
4216 !C Varibles set size of box
4217 with_theta_constr = index(controlcard,"WITH_THETA_CONSTR").gt.0
4218 protein=index(controlcard,"PROTEIN").gt.0
4219 ions=index(controlcard,"IONS").gt.0
4220 nucleic=index(controlcard,"NUCLEIC").gt.0
4221 write (iout,*) "with_theta_constr ",with_theta_constr
4222 AFMlog=(index(controlcard,'AFM'))
4223 selfguide=(index(controlcard,'SELFGUIDE'))
4224 print *,'AFMlog',AFMlog,selfguide,"KUPA"
4225 call readi(controlcard,'GENCONSTR',genconstr,0)
4226 call reada(controlcard,'BOXX',boxxsize,100.0d0)
4227 call reada(controlcard,'BOXY',boxysize,100.0d0)
4228 call reada(controlcard,'BOXZ',boxzsize,100.0d0)
4229 call readi(controlcard,'TUBEMOD',tubemode,0)
4230 print *,"SCELE",scelemode
4231 call readi(controlcard,"SCELEMODE",scelemode,0)
4232 print *,"SCELE",scelemode
4234 ! elemode = 0 is orignal UNRES electrostatics
4235 ! elemode = 1 is "Momo" potentials in progress
4236 ! elemode = 2 is in development EVALD
4237 write (iout,*) TUBEmode,"TUBEMODE"
4238 if (TUBEmode.gt.0) then
4239 call reada(controlcard,"XTUBE",tubecenter(1),0.0d0)
4240 call reada(controlcard,"YTUBE",tubecenter(2),0.0d0)
4241 call reada(controlcard,"ZTUBE",tubecenter(3),0.0d0)
4242 call reada(controlcard,"RTUBE",tubeR0,0.0d0)
4243 call reada(controlcard,"TUBETOP",bordtubetop,boxzsize)
4244 call reada(controlcard,"TUBEBOT",bordtubebot,0.0d0)
4245 call reada(controlcard,"TUBEBUF",tubebufthick,1.0d0)
4246 buftubebot=bordtubebot+tubebufthick
4247 buftubetop=bordtubetop-tubebufthick
4250 ! CUTOFFF ON ELECTROSTATICS
4251 call reada(controlcard,"R_CUT_ELE",r_cut_ele,15.0d0)
4252 call reada(controlcard,"LAMBDA_ELE",rlamb_ele,0.3d0)
4253 write(iout,*) "R_CUT_ELE=",r_cut_ele
4254 ! Lipidic parameters
4255 call reada(controlcard,"LIPTHICK",lipthick,0.0d0)
4256 call reada(controlcard,"LIPAQBUF",lipbufthick,0.0d0)
4257 if (lipthick.gt.0.0d0) then
4258 bordliptop=(boxzsize+lipthick)/2.0
4259 bordlipbot=bordliptop-lipthick
4260 if ((bordliptop.gt.boxzsize).or.(bordlipbot.lt.0.0)) &
4261 write(iout,*) "WARNING WRONG SIZE OF LIPIDIC PHASE"
4262 buflipbot=bordlipbot+lipbufthick
4263 bufliptop=bordliptop-lipbufthick
4264 if ((lipbufthick*2.0d0).gt.lipthick) &
4265 write(iout,*) "WARNING WRONG SIZE OF LIP AQ BUF"
4266 endif !lipthick.gt.0
4267 write(iout,*) "bordliptop=",bordliptop
4268 write(iout,*) "bordlipbot=",bordlipbot
4269 write(iout,*) "bufliptop=",bufliptop
4270 write(iout,*) "buflipbot=",buflipbot
4271 write (iout,*) "SHIELD MODE",shield_mode
4273 !C-------------------------
4274 minim=(index(controlcard,'MINIMIZE').gt.0)
4275 dccart=(index(controlcard,'CART').gt.0)
4276 overlapsc=(index(controlcard,'OVERLAP').gt.0)
4277 overlapsc=.not.overlapsc
4278 searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
4279 searchsc=.not.searchsc
4280 sideadd=(index(controlcard,'SIDEADD').gt.0)
4281 energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
4282 outpdb=(index(controlcard,'PDBOUT').gt.0)
4283 outmol2=(index(controlcard,'MOL2OUT').gt.0)
4284 pdbref=(index(controlcard,'PDBREF').gt.0)
4285 refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
4286 indpdb=index(controlcard,'PDBSTART')
4287 extconf=(index(controlcard,'EXTCONF').gt.0)
4288 call readi(controlcard,'IPRINT',iprint,0)
4289 call readi(controlcard,'MAXGEN',maxgen,10000)
4290 call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
4291 call readi(controlcard,"KDIAG",kdiag,0)
4292 call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
4293 if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) &
4294 write (iout,*) "RESCALE_MODE",rescale_mode
4295 split_ene=index(controlcard,'SPLIT_ENE').gt.0
4296 if (index(controlcard,'REGULAR').gt.0.0D0) then
4297 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4301 if (index(controlcard,'CHECKGRAD').gt.0) then
4303 if (index(controlcard,'CART').gt.0) then
4305 elseif (index(controlcard,'CARINT').gt.0) then
4310 elseif (index(controlcard,'THREAD').gt.0) then
4312 call readi(controlcard,'THREAD',nthread,0)
4313 if (nthread.gt.0) then
4314 call reada(controlcard,'WEIDIS',weidis,0.1D0)
4317 write (iout,'(a)')'A number has to follow the THREAD keyword.'
4318 stop 'Error termination in Read_Control.'
4320 else if (index(controlcard,'MCMA').gt.0) then
4322 else if (index(controlcard,'MCEE').gt.0) then
4324 else if (index(controlcard,'MULTCONF').gt.0) then
4326 else if (index(controlcard,'MAP').gt.0) then
4328 call readi(controlcard,'MAP',nmap,0)
4329 else if (index(controlcard,'CSA').gt.0) then
4331 !rc else if (index(controlcard,'ZSCORE').gt.0) then
4333 !rc ZSCORE is rm from UNRES, modecalc=9 is available
4336 !fcm else if (index(controlcard,'MCMF').gt.0) then
4338 else if (index(controlcard,'SOFTREG').gt.0) then
4340 else if (index(controlcard,'CHECK_BOND').gt.0) then
4342 else if (index(controlcard,'TEST').gt.0) then
4344 else if (index(controlcard,'MD').gt.0) then
4346 else if (index(controlcard,'RE ').gt.0) then
4350 lmuca=index(controlcard,'MUCA').gt.0
4351 call readi(controlcard,'MUCADYN',mucadyn,0)
4352 call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
4353 if (lmuca .and. (me.eq.king .or. .not.out1file )) &
4355 write (iout,*) 'MUCADYN=',mucadyn
4356 write (iout,*) 'MUCASMOOTH=',muca_smooth
4359 iscode=index(controlcard,'ONE_LETTER')
4360 indphi=index(controlcard,'PHI')
4361 indback=index(controlcard,'BACK')
4362 iranconf=index(controlcard,'RAND_CONF')
4363 i2ndstr=index(controlcard,'USE_SEC_PRED')
4364 gradout=index(controlcard,'GRADOUT').gt.0
4365 gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
4366 call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
4367 if (me.eq.king .or. .not.out1file ) &
4368 write (iout,*) "DISTCHAINMAX",distchainmax
4370 if(me.eq.king.or..not.out1file) &
4371 write (iout,'(2a)') diagmeth(kdiag),&
4372 ' routine used to diagonalize matrices.'
4373 if (shield_mode.gt.0) then
4374 pi=4.0D0*datan(1.0D0)
4375 !C VSolvSphere the volume of solving sphere
4377 !C rpp(1,1) is the energy r0 for peptide group contact and will be used for it
4378 !C there will be no distinction between proline peptide group and normal peptide
4379 !C group in case of shielding parameters
4380 VSolvSphere=4.0/3.0*pi*(4.50d0)**3
4381 VSolvSphere_div=VSolvSphere-4.0/3.0*pi*(4.50/2.0)**3
4382 write (iout,*) VSolvSphere,VSolvSphere_div
4383 !C long axis of side chain
4385 ! long_r_sidechain(i)=vbldsc0(1,i)
4386 ! short_r_sidechain(i)=sigma0(i)
4387 ! write(iout,*) "parame for long and short axis",i,vbldsc0(1,i),&
4393 end subroutine read_control
4394 !-----------------------------------------------------------------------------
4395 subroutine read_REMDpar
4397 ! Read REMD settings
4404 use control_data, only:out1file
4405 ! implicit real*8 (a-h,o-z)
4406 ! include 'DIMENSIONS'
4407 ! include 'COMMON.IOUNITS'
4408 ! include 'COMMON.TIME1'
4409 ! include 'COMMON.MD'
4412 !el include 'COMMON.LANGEVIN'
4414 !el include 'COMMON.LANGEVIN.lang0'
4416 ! include 'COMMON.INTERACT'
4417 ! include 'COMMON.NAMES'
4418 ! include 'COMMON.GEO'
4419 ! include 'COMMON.REMD'
4420 ! include 'COMMON.CONTROL'
4421 ! include 'COMMON.SETUP'
4422 ! character(len=80) :: ucase
4423 character(len=320) :: controlcard
4424 character(len=3200) :: controlcard1
4425 integer :: iremd_m_total
4428 ! real(kind=8) :: var,ene
4430 if(me.eq.king.or..not.out1file) &
4431 write (iout,*) "REMD setup"
4433 call card_concat(controlcard,.true.)
4434 call readi(controlcard,"NREP",nrep,3)
4435 call readi(controlcard,"NSTEX",nstex,1000)
4436 call reada(controlcard,"RETMIN",retmin,10.0d0)
4437 call reada(controlcard,"RETMAX",retmax,1000.0d0)
4438 mremdsync=(index(controlcard,'SYNC').gt.0)
4439 call readi(controlcard,"NSYN",i_sync_step,100)
4440 restart1file=(index(controlcard,'REST1FILE').gt.0)
4441 traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
4442 call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
4443 if(max_cache_traj_use.gt.max_cache_traj) &
4444 max_cache_traj_use=max_cache_traj
4445 if(me.eq.king.or..not.out1file) then
4446 !d if (traj1file) then
4447 !rc caching is in testing - NTWX is not ignored
4448 !d write (iout,*) "NTWX value is ignored"
4449 !d write (iout,*) " trajectory is stored to one file by master"
4450 !d write (iout,*) " before exchange at NSTEX intervals"
4452 write (iout,*) "NREP= ",nrep
4453 write (iout,*) "NSTEX= ",nstex
4454 write (iout,*) "SYNC= ",mremdsync
4455 write (iout,*) "NSYN= ",i_sync_step
4456 write (iout,*) "TRAJCACHE= ",max_cache_traj_use
4459 allocate(remd_t(nrep),remd_m(nrep)) !(maxprocs)
4460 if (index(controlcard,'TLIST').gt.0) then
4462 call card_concat(controlcard1,.true.)
4463 read(controlcard1,*) (remd_t(i),i=1,nrep)
4464 if(me.eq.king.or..not.out1file) &
4465 write (iout,*)'tlist',(remd_t(i),i=1,nrep)
4468 if (index(controlcard,'MLIST').gt.0) then
4470 call card_concat(controlcard1,.true.)
4471 read(controlcard1,*) (remd_m(i),i=1,nrep)
4472 if(me.eq.king.or..not.out1file) then
4473 write (iout,*)'mlist',(remd_m(i),i=1,nrep)
4476 iremd_m_total=iremd_m_total+remd_m(i)
4478 write (iout,*) 'Total number of replicas ',iremd_m_total
4481 if(me.eq.king.or..not.out1file) &
4482 write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
4484 end subroutine read_REMDpar
4485 !-----------------------------------------------------------------------------
4486 subroutine read_MDpar
4490 use control_data, only: r_cut,rlamb,out1file
4492 use geometry_data, only: pi
4494 ! implicit real*8 (a-h,o-z)
4495 ! include 'DIMENSIONS'
4496 ! include 'COMMON.IOUNITS'
4497 ! include 'COMMON.TIME1'
4498 ! include 'COMMON.MD'
4501 !el include 'COMMON.LANGEVIN'
4503 !el include 'COMMON.LANGEVIN.lang0'
4505 ! include 'COMMON.INTERACT'
4506 ! include 'COMMON.NAMES'
4507 ! include 'COMMON.GEO'
4508 ! include 'COMMON.SETUP'
4509 ! include 'COMMON.CONTROL'
4510 ! include 'COMMON.SPLITELE'
4511 ! character(len=80) :: ucase
4512 character(len=320) :: controlcard
4517 call card_concat(controlcard,.true.)
4518 call readi(controlcard,"NSTEP",n_timestep,1000000)
4519 call readi(controlcard,"NTWE",ntwe,100)
4520 call readi(controlcard,"NTWX",ntwx,1000)
4521 call reada(controlcard,"DT",d_time,1.0d-1)
4522 call reada(controlcard,"DVMAX",dvmax,2.0d1)
4523 call reada(controlcard,"DAMAX",damax,1.0d1)
4524 call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
4525 call readi(controlcard,"LANG",lang,0)
4526 RESPA = index(controlcard,"RESPA") .gt. 0
4527 call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
4528 ntime_split0=ntime_split
4529 call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
4530 ntime_split0=ntime_split
4531 call reada(controlcard,"R_CUT",r_cut,2.0d0)
4532 call reada(controlcard,"LAMBDA",rlamb,0.3d0)
4533 rest = index(controlcard,"REST").gt.0
4534 tbf = index(controlcard,"TBF").gt.0
4535 usampl = index(controlcard,"USAMPL").gt.0
4536 mdpdb = index(controlcard,"MDPDB").gt.0
4537 call reada(controlcard,"T_BATH",t_bath,300.0d0)
4538 call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
4539 call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
4540 call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
4541 if (count_reset_moment.eq.0) count_reset_moment=1000000000
4542 call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
4543 reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
4544 reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
4545 if (count_reset_vel.eq.0) count_reset_vel=1000000000
4546 large = index(controlcard,"LARGE").gt.0
4547 print_compon = index(controlcard,"PRINT_COMPON").gt.0
4548 rattle = index(controlcard,"RATTLE").gt.0
4549 preminim=(index(controlcard,'PREMINIM').gt.0)
4550 write (iout,*) "PREMINIM ",preminim
4551 dccart=(index(controlcard,'CART').gt.0)
4552 if (preminim) call read_minim
4553 ! if performing umbrella sampling, fragments constrained are read from the fragment file
4559 if(me.eq.king.or..not.out1file) then
4561 write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
4563 write (iout,'(a)') "The units are:"
4564 write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
4565 write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",&
4566 " acceleration: angstrom/(48.9 fs)**2"
4567 write (iout,'(a)') "energy: kcal/mol, temperature: K"
4569 write (iout,'(a60,i10)') "Number of time steps:",n_timestep
4570 write (iout,'(a60,f10.5,a)') &
4571 "Initial time step of numerical integration:",d_time,&
4573 write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
4575 write (iout,'(2a,i4,a)') &
4576 "A-MTS algorithm used; initial time step for fast-varying",&
4577 " short-range forces split into",ntime_split," steps."
4578 write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",&
4579 r_cut," lambda",rlamb
4581 write (iout,'(2a,f10.5)') &
4582 "Maximum acceleration threshold to reduce the time step",&
4583 "/increase split number:",damax
4584 write (iout,'(2a,f10.5)') &
4585 "Maximum predicted energy drift to reduce the timestep",&
4586 "/increase split number:",edriftmax
4587 write (iout,'(a60,f10.5)') &
4588 "Maximum velocity threshold to reduce velocities:",dvmax
4589 write (iout,'(a60,i10)') "Frequency of property output:",ntwe
4590 write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
4591 if (rattle) write (iout,'(a60)') &
4592 "Rattle algorithm used to constrain the virtual bonds"
4596 call reada(controlcard,"ETAWAT",etawat,0.8904d0)
4597 call reada(controlcard,"RWAT",rwat,1.4d0)
4598 call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
4599 surfarea=index(controlcard,"SURFAREA").gt.0
4600 call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
4601 if(me.eq.king.or..not.out1file)then
4602 write (iout,'(/a,$)') "Langevin dynamics calculation"
4604 write (iout,'(a/)') &
4605 " with direct integration of Langevin equations"
4606 else if (lang.eq.2) then
4607 write (iout,'(a/)') " with TINKER stochasic MD integrator"
4608 else if (lang.eq.3) then
4609 write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
4610 else if (lang.eq.4) then
4611 write (iout,'(a/)') " in overdamped mode"
4613 write (iout,'(//a,i5)') &
4614 "=========== ERROR: Unknown Langevin dynamics mode:",lang
4617 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4618 write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
4619 write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
4620 write (iout,'(a60,f10.5)') &
4621 "Scaling factor of the friction forces:",scal_fric
4622 if (surfarea) write (iout,'(2a,i10,a)') &
4623 "Friction coefficients will be scaled by solvent-accessible",&
4624 " surface area every",reset_fricmat," steps."
4626 ! Calculate friction coefficients and bounds of stochastic forces
4627 eta=6*pi*cPoise*etawat
4628 if(me.eq.king.or..not.out1file) &
4629 write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:",&
4632 do j=1,5 !types of molecules
4633 gamp(j)=scal_fric*(pstok(j)+rwat)*eta
4634 stdfp(j)=dsqrt(2*Rb*t_bath/d_time)
4636 allocate(gamsc(ntyp1,5),stdfsc(ntyp1,5)) !(ntyp1)
4637 do j=1,5 !types of molecules
4639 gamsc(i,j)=scal_fric*(restok(i,j)+rwat)*eta
4640 stdfsc(i,j)=dsqrt(2*Rb*t_bath/d_time)
4644 if(me.eq.king.or..not.out1file)then
4645 write (iout,'(/2a/)') &
4646 "Radii of site types and friction coefficients and std's of",&
4647 " stochastic forces of fully exposed sites"
4648 write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp(1),stdfp*dsqrt(gamp(1))
4650 write (iout,'(a5,f5.2,2f10.5)') restyp(i,1),restok(i,1),&
4651 gamsc(i,1),stdfsc(i,1)*dsqrt(gamsc(i,1))
4655 if(me.eq.king.or..not.out1file)then
4656 write (iout,'(a)') "Berendsen bath calculation"
4657 write (iout,'(a60,f10.5)') "Temperature:",t_bath
4658 write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
4660 write (iout,'(a,i10,a)') "Momenta will be reset at zero every",&
4661 count_reset_moment," steps"
4663 write (iout,'(a,i10,a)') &
4664 "Velocities will be reset at random every",count_reset_vel,&
4668 if(me.eq.king.or..not.out1file) &
4669 write (iout,'(a31)') "Microcanonical mode calculation"
4671 if(me.eq.king.or..not.out1file)then
4672 if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
4674 write(iout,*) "MD running with constraints."
4675 write(iout,*) "Equilibration time ", eq_time, " mtus."
4676 write(iout,*) "Constraining ", nfrag," fragments."
4677 write(iout,*) "Length of each fragment, weight and q0:"
4679 write (iout,*) "Set of restraints #",iset
4681 write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),&
4682 ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
4684 write(iout,*) "constraints between ", npair, "fragments."
4685 write(iout,*) "constraint pairs, weights and q0:"
4687 write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),&
4688 ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
4690 write(iout,*) "angle constraints within ", nfrag_back,&
4691 "backbone fragments."
4692 write(iout,*) "fragment, weights:"
4694 write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),&
4695 ifrag_back(2,i,iset),wfrag_back(1,i,iset),&
4696 wfrag_back(2,i,iset),wfrag_back(3,i,iset)
4699 iset=mod(kolor,nset)+1
4702 if(me.eq.king.or..not.out1file) &
4703 write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
4705 end subroutine read_MDpar
4706 !-----------------------------------------------------------------------------
4710 ! implicit real*8 (a-h,o-z)
4711 ! include 'DIMENSIONS'
4712 ! include 'COMMON.MAP'
4713 ! include 'COMMON.IOUNITS'
4714 character(len=3) :: angid(4) = (/'THE','PHI','ALP','OME'/)
4715 character(len=80) :: mapcard !,ucase
4718 ! real(kind=8) :: var,ene
4721 read (inp,'(a)') mapcard
4722 mapcard=ucase(mapcard)
4723 if (index(mapcard,'PHI').gt.0) then
4725 else if (index(mapcard,'THE').gt.0) then
4727 else if (index(mapcard,'ALP').gt.0) then
4729 else if (index(mapcard,'OME').gt.0) then
4732 write(iout,'(a)')'Error - illegal variable spec in MAP card.'
4733 stop 'Error - illegal variable spec in MAP card.'
4735 call readi (mapcard,'RES1',res1(imap),0)
4736 call readi (mapcard,'RES2',res2(imap),0)
4737 if (res1(imap).eq.0) then
4738 res1(imap)=res2(imap)
4739 else if (res2(imap).eq.0) then
4740 res2(imap)=res1(imap)
4742 if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
4743 write (iout,'(a)') &
4744 'Error - illegal definition of variable group in MAP.'
4745 stop 'Error - illegal definition of variable group in MAP.'
4747 call reada(mapcard,'FROM',ang_from(imap),0.0D0)
4748 call reada(mapcard,'TO',ang_to(imap),0.0D0)
4749 call readi(mapcard,'NSTEP',nstep(imap),0)
4750 if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
4751 write (iout,'(a)') &
4752 'Illegal boundary and/or step size specification in MAP.'
4753 stop 'Illegal boundary and/or step size specification in MAP.'
4757 end subroutine map_read
4758 !-----------------------------------------------------------------------------
4761 use control_data, only: vdisulf
4763 ! implicit real*8 (a-h,o-z)
4764 ! include 'DIMENSIONS'
4765 ! include 'COMMON.IOUNITS'
4766 ! include 'COMMON.GEO'
4767 ! include 'COMMON.CSA'
4768 ! include 'COMMON.BANK'
4769 ! include 'COMMON.CONTROL'
4770 ! character(len=80) :: ucase
4771 character(len=620) :: mcmcard
4773 ! integer :: ntf,ik,iw_pdb
4774 ! real(kind=8) :: var,ene
4776 call card_concat(mcmcard,.true.)
4778 call readi(mcmcard,'NCONF',nconf,50)
4779 call readi(mcmcard,'NADD',nadd,0)
4780 call readi(mcmcard,'JSTART',jstart,1)
4781 call readi(mcmcard,'JEND',jend,1)
4782 call readi(mcmcard,'NSTMAX',nstmax,500000)
4783 call readi(mcmcard,'N0',n0,1)
4784 call readi(mcmcard,'N1',n1,6)
4785 call readi(mcmcard,'N2',n2,4)
4786 call readi(mcmcard,'N3',n3,0)
4787 call readi(mcmcard,'N4',n4,0)
4788 call readi(mcmcard,'N5',n5,0)
4789 call readi(mcmcard,'N6',n6,10)
4790 call readi(mcmcard,'N7',n7,0)
4791 call readi(mcmcard,'N8',n8,0)
4792 call readi(mcmcard,'N9',n9,0)
4793 call readi(mcmcard,'N14',n14,0)
4794 call readi(mcmcard,'N15',n15,0)
4795 call readi(mcmcard,'N16',n16,0)
4796 call readi(mcmcard,'N17',n17,0)
4797 call readi(mcmcard,'N18',n18,0)
4799 vdisulf=(index(mcmcard,'DYNSS').gt.0)
4801 call readi(mcmcard,'NDIFF',ndiff,2)
4802 call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
4803 call readi(mcmcard,'IS1',is1,1)
4804 call readi(mcmcard,'IS2',is2,8)
4805 call readi(mcmcard,'NRAN0',nran0,4)
4806 call readi(mcmcard,'NRAN1',nran1,2)
4807 call readi(mcmcard,'IRR',irr,1)
4808 call readi(mcmcard,'NSEED',nseed,20)
4809 call readi(mcmcard,'NTOTAL',ntotal,10000)
4810 call reada(mcmcard,'CUT1',cut1,2.0d0)
4811 call reada(mcmcard,'CUT2',cut2,5.0d0)
4812 call reada(mcmcard,'ESTOP',estop,-3000.0d0)
4813 call readi(mcmcard,'ICMAX',icmax,3)
4814 call readi(mcmcard,'IRESTART',irestart,0)
4815 !!bankt call readi(mcmcard,'NBANKTM',ntbankm,0)
4818 call reada(mcmcard,'DELE',dele,20.0d0)
4819 call reada(mcmcard,'DIFCUT',difcut,720.0d0)
4820 call readi(mcmcard,'IREF',iref,0)
4821 call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
4822 call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
4823 call readi(mcmcard,'NCONF_IN',nconf_in,0)
4824 call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
4825 write (iout,*) "NCONF_IN",nconf_in
4827 end subroutine csaread
4828 !-----------------------------------------------------------------------------
4832 use control_data, only: MaxMoveType
4835 ! implicit real*8 (a-h,o-z)
4836 ! include 'DIMENSIONS'
4837 ! include 'COMMON.MCM'
4838 ! include 'COMMON.MCE'
4839 ! include 'COMMON.IOUNITS'
4840 ! character(len=80) :: ucase
4841 character(len=320) :: mcmcard
4844 ! real(kind=8) :: var,ene
4846 call card_concat(mcmcard,.true.)
4847 call readi(mcmcard,'MAXACC',maxacc,100)
4848 call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
4849 call readi(mcmcard,'MAXTRIAL',maxtrial,100)
4850 call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
4851 call readi(mcmcard,'MAXREPM',maxrepm,200)
4852 call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
4853 call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
4854 call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
4855 call reada(mcmcard,'E_UP',e_up,5.0D0)
4856 call reada(mcmcard,'DELTE',delte,0.1D0)
4857 call readi(mcmcard,'NSWEEP',nsweep,5)
4858 call readi(mcmcard,'NSTEPH',nsteph,0)
4859 call readi(mcmcard,'NSTEPC',nstepc,0)
4860 call reada(mcmcard,'TMIN',tmin,298.0D0)
4861 call reada(mcmcard,'TMAX',tmax,298.0D0)
4862 call readi(mcmcard,'NWINDOW',nwindow,0)
4863 call readi(mcmcard,'PRINT_MC',print_mc,0)
4864 print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
4865 print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
4866 ent_read=(index(mcmcard,'ENT_READ').gt.0)
4867 call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
4868 call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
4869 call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
4870 call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
4871 call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
4872 if (nwindow.gt.0) then
4873 allocate(winstart(nwindow)) !!el (maxres)
4874 allocate(winend(nwindow)) !!el
4875 allocate(winlen(nwindow)) !!el
4876 read (inp,*) (winstart(i),winend(i),i=1,nwindow)
4878 winlen(i)=winend(i)-winstart(i)+1
4881 if (tmax.lt.tmin) tmax=tmin
4882 if (tmax.eq.tmin) then
4886 if (nstepc.gt.0 .and. nsteph.gt.0) then
4887 tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
4888 tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
4890 allocate(sumpro_type(0:MaxMoveType)) !(0:MaxMoveType)
4891 ! Probabilities of different move types
4892 sumpro_type(0)=0.0D0
4893 call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
4894 call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
4895 sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
4896 call reada(mcmcard,'THETA' ,sumpro_type(3),0.0d0)
4897 sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
4898 call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
4899 sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
4901 print *,'i',i,' sumprotype',sumpro_type(i)
4902 sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
4903 print *,'i',i,' sumprotype',sumpro_type(i)
4906 end subroutine mcmread
4907 !-----------------------------------------------------------------------------
4908 subroutine read_minim
4911 ! implicit real*8 (a-h,o-z)
4912 ! include 'DIMENSIONS'
4913 ! include 'COMMON.MINIM'
4914 ! include 'COMMON.IOUNITS'
4915 ! character(len=80) :: ucase
4916 character(len=320) :: minimcard
4918 ! integer :: ntf,ik,iw_pdb
4919 ! real(kind=8) :: var,ene
4921 call card_concat(minimcard,.true.)
4922 call readi(minimcard,'MAXMIN',maxmin,2000)
4923 call readi(minimcard,'MAXFUN',maxfun,5000)
4924 call readi(minimcard,'MINMIN',minmin,maxmin)
4925 call readi(minimcard,'MINFUN',minfun,maxmin)
4926 call reada(minimcard,'TOLF',tolf,1.0D-2)
4927 call reada(minimcard,'RTOLF',rtolf,1.0D-4)
4928 print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
4929 print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
4930 print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
4931 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
4932 'Options in energy minimization:'
4933 write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)') &
4934 'MaxMin:',MaxMin,' MaxFun:',MaxFun,&
4935 'MinMin:',MinMin,' MinFun:',MinFun,&
4936 ' TolF:',TolF,' RTolF:',RTolF
4938 end subroutine read_minim
4939 !-----------------------------------------------------------------------------
4940 subroutine openunits
4942 use MD_data, only: usampl
4945 use control_data, only:out1file
4946 use control, only: getenv_loc
4947 ! implicit real*8 (a-h,o-z)
4948 ! include 'DIMENSIONS'
4951 character(len=16) :: form,nodename
4952 integer :: nodelen,ierror,npos
4954 ! include 'COMMON.SETUP'
4955 ! include 'COMMON.IOUNITS'
4956 ! include 'COMMON.MD'
4957 ! include 'COMMON.CONTROL'
4958 integer :: lenpre,lenpot,lentmp !,ilen
4960 character(len=3) :: out1file_text !,ucase
4961 character(len=3) :: ll
4964 ! integer :: ntf,ik,iw_pdb
4965 ! real(kind=8) :: var,ene
4967 ! print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
4968 call getenv_loc("PREFIX",prefix)
4970 call getenv_loc("POT",pot)
4971 call getenv_loc("DIRTMP",tmpdir)
4972 call getenv_loc("CURDIR",curdir)
4973 call getenv_loc("OUT1FILE",out1file_text)
4974 ! print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
4975 out1file_text=ucase(out1file_text)
4976 if (out1file_text(1:1).eq."Y") then
4979 out1file=fg_rank.gt.0
4984 if (lentmp.gt.0) then
4985 write (*,'(80(1h!))')
4986 write (*,'(a,19x,a,19x,a)') "!"," A T T E N T I O N ","!"
4987 write (*,'(80(1h!))')
4988 write (*,*)"All output files will be on node /tmp directory."
4990 call MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
4991 if (me.eq.king) then
4992 write (*,*) "The master node is ",nodename
4993 else if (fg_rank.eq.0) then
4994 write (*,*) "I am the CG slave node ",nodename
4996 write (*,*) "I am the FG slave node ",nodename
4999 PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
5000 lenpre = lentmp+lenpre+1
5002 entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
5003 ! Get the names and open the input files
5004 #if defined(WINIFL) || defined(WINPGI)
5005 open(1,file=pref_orig(:ilen(pref_orig))// &
5006 '.inp',status='old',readonly,shared)
5007 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5008 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5009 ! Get parameter filenames and open the parameter files.
5010 call getenv_loc('BONDPAR',bondname)
5011 open (ibond,file=bondname,status='old',readonly,shared)
5012 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5013 open (ibond_nucl,file=bondname_nucl,status='old',readonly,shared)
5014 call getenv_loc('THETPAR',thetname)
5015 open (ithep,file=thetname,status='old',readonly,shared)
5016 call getenv_loc('ROTPAR',rotname)
5017 open (irotam,file=rotname,status='old',readonly,shared)
5018 call getenv_loc('TORPAR',torname)
5019 open (itorp,file=torname,status='old',readonly,shared)
5020 call getenv_loc('TORDPAR',tordname)
5021 open (itordp,file=tordname,status='old',readonly,shared)
5022 call getenv_loc('FOURIER',fouriername)
5023 open (ifourier,file=fouriername,status='old',readonly,shared)
5024 call getenv_loc('ELEPAR',elename)
5025 open (ielep,file=elename,status='old',readonly,shared)
5026 call getenv_loc('SIDEPAR',sidename)
5027 open (isidep,file=sidename,status='old',readonly,shared)
5029 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5030 open (ithep_nucl,file=thetname_nucl,status='old',readonly,shared)
5031 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5032 open (irotam_nucl,file=rotname_nucl,status='old',readonly,shared)
5033 call getenv_loc('TORPAR_NUCL',torname_nucl)
5034 open (itorp_nucl,file=torname_nucl,status='old',readonly,shared)
5035 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5036 open (itordp_nucl,file=tordname_nucl,status='old',readonly,shared)
5037 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5038 open (isidep_nucl,file=sidename_nucl,status='old',readonly,shared)
5041 #elif (defined CRAY) || (defined AIX)
5042 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5044 ! print *,"Processor",myrank," opened file 1"
5045 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5046 ! print *,"Processor",myrank," opened file 9"
5047 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5048 ! Get parameter filenames and open the parameter files.
5049 call getenv_loc('BONDPAR',bondname)
5050 open (ibond,file=bondname,status='old',action='read')
5051 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5052 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5054 ! print *,"Processor",myrank," opened file IBOND"
5055 call getenv_loc('THETPAR',thetname)
5056 open (ithep,file=thetname,status='old',action='read')
5057 ! print *,"Processor",myrank," opened file ITHEP"
5058 call getenv_loc('ROTPAR',rotname)
5059 open (irotam,file=rotname,status='old',action='read')
5060 ! print *,"Processor",myrank," opened file IROTAM"
5061 call getenv_loc('TORPAR',torname)
5062 open (itorp,file=torname,status='old',action='read')
5063 ! print *,"Processor",myrank," opened file ITORP"
5064 call getenv_loc('TORDPAR',tordname)
5065 open (itordp,file=tordname,status='old',action='read')
5066 ! print *,"Processor",myrank," opened file ITORDP"
5067 call getenv_loc('SCCORPAR',sccorname)
5068 open (isccor,file=sccorname,status='old',action='read')
5069 ! print *,"Processor",myrank," opened file ISCCOR"
5070 call getenv_loc('FOURIER',fouriername)
5071 open (ifourier,file=fouriername,status='old',action='read')
5072 ! print *,"Processor",myrank," opened file IFOURIER"
5073 call getenv_loc('ELEPAR',elename)
5074 open (ielep,file=elename,status='old',action='read')
5075 ! print *,"Processor",myrank," opened file IELEP"
5076 call getenv_loc('SIDEPAR',sidename)
5077 open (isidep,file=sidename,status='old',action='read')
5079 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5080 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5081 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5082 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5083 call getenv_loc('TORPAR_NUCL',torname_nucl)
5084 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5085 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5086 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5087 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5088 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5090 call getenv_loc('LIPTRANPAR',liptranname)
5091 open (iliptranpar,file=liptranname,status='old',action='read')
5092 call getenv_loc('TUBEPAR',tubename)
5093 open (itube,file=tubename,status='old',action='read')
5094 call getenv_loc('IONPAR',ionname)
5095 open (iion,file=ionname,status='old',action='read')
5097 ! print *,"Processor",myrank," opened file ISIDEP"
5098 ! print *,"Processor",myrank," opened parameter files"
5100 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
5101 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5102 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5103 ! Get parameter filenames and open the parameter files.
5104 call getenv_loc('BONDPAR',bondname)
5105 open (ibond,file=bondname,status='old')
5106 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5107 open (ibond_nucl,file=bondname_nucl,status='old')
5109 call getenv_loc('THETPAR',thetname)
5110 open (ithep,file=thetname,status='old')
5111 call getenv_loc('ROTPAR',rotname)
5112 open (irotam,file=rotname,status='old')
5113 call getenv_loc('TORPAR',torname)
5114 open (itorp,file=torname,status='old')
5115 call getenv_loc('TORDPAR',tordname)
5116 open (itordp,file=tordname,status='old')
5117 call getenv_loc('SCCORPAR',sccorname)
5118 open (isccor,file=sccorname,status='old')
5119 call getenv_loc('FOURIER',fouriername)
5120 open (ifourier,file=fouriername,status='old')
5121 call getenv_loc('ELEPAR',elename)
5122 open (ielep,file=elename,status='old')
5123 call getenv_loc('SIDEPAR',sidename)
5124 open (isidep,file=sidename,status='old')
5126 open (ithep_nucl,file=thetname_nucl,status='old')
5127 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5128 open (irotam_nucl,file=rotname_nucl,status='old')
5129 call getenv_loc('TORPAR_NUCL',torname_nucl)
5130 open (itorp_nucl,file=torname_nucl,status='old')
5131 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5132 open (itordp_nucl,file=tordname_nucl,status='old')
5133 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5134 open (isidep_nucl,file=sidename_nucl,status='old')
5136 call getenv_loc('LIPTRANPAR',liptranname)
5137 open (iliptranpar,file=liptranname,status='old')
5138 call getenv_loc('TUBEPAR',tubename)
5139 open (itube,file=tubename,status='old')
5140 call getenv_loc('IONPAR',ionname)
5141 open (iion,file=ionname,status='old')
5143 open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',&
5145 open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
5146 ! open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
5147 ! Get parameter filenames and open the parameter files.
5148 call getenv_loc('BONDPAR',bondname)
5149 open (ibond,file=bondname,status='old',action='read')
5150 call getenv_loc('BONDPAR_NUCL',bondname_nucl)
5151 open (ibond_nucl,file=bondname_nucl,status='old',action='read')
5152 call getenv_loc('THETPAR',thetname)
5153 open (ithep,file=thetname,status='old',action='read')
5154 call getenv_loc('ROTPAR',rotname)
5155 open (irotam,file=rotname,status='old',action='read')
5156 call getenv_loc('TORPAR',torname)
5157 open (itorp,file=torname,status='old',action='read')
5158 call getenv_loc('TORDPAR',tordname)
5159 open (itordp,file=tordname,status='old',action='read')
5160 call getenv_loc('SCCORPAR',sccorname)
5161 open (isccor,file=sccorname,status='old',action='read')
5163 call getenv_loc('THETPARPDB',thetname_pdb)
5164 print *,"thetname_pdb ",thetname_pdb
5165 open (ithep_pdb,file=thetname_pdb,status='old',action='read')
5166 print *,ithep_pdb," opened"
5168 call getenv_loc('FOURIER',fouriername)
5169 open (ifourier,file=fouriername,status='old',readonly)
5170 call getenv_loc('ELEPAR',elename)
5171 open (ielep,file=elename,status='old',readonly)
5172 call getenv_loc('SIDEPAR',sidename)
5173 open (isidep,file=sidename,status='old',readonly)
5175 call getenv_loc('THETPAR_NUCL',thetname_nucl)
5176 open (ithep_nucl,file=thetname_nucl,status='old',action='read')
5177 call getenv_loc('ROTPAR_NUCL',rotname_nucl)
5178 open (irotam_nucl,file=rotname_nucl,status='old',action='read')
5179 call getenv_loc('TORPAR_NUCL',torname_nucl)
5180 open (itorp_nucl,file=torname_nucl,status='old',action='read')
5181 call getenv_loc('TORDPAR_NUCL',tordname_nucl)
5182 open (itordp_nucl,file=tordname_nucl,status='old',action='read')
5183 call getenv_loc('SIDEPAR_NUCL',sidename_nucl)
5184 open (isidep_nucl,file=sidename_nucl,status='old',action='read')
5185 call getenv_loc('SIDEPAR_SCBASE',sidename_scbase)
5186 open (isidep_scbase,file=sidename_scbase,status='old',action='read')
5187 call getenv_loc('PEPPAR_PEPBASE',pepname_pepbase)
5188 open (isidep_pepbase,file=pepname_pepbase,status='old',action='read')
5189 call getenv_loc('SCPAR_PHOSPH',pepname_scpho)
5190 open (isidep_scpho,file=pepname_scpho,status='old',action='read')
5191 call getenv_loc('PEPPAR_PHOSPH',pepname_peppho)
5192 open (isidep_peppho,file=pepname_peppho,status='old',action='read')
5195 call getenv_loc('LIPTRANPAR',liptranname)
5196 open (iliptranpar,file=liptranname,status='old',action='read')
5197 call getenv_loc('TUBEPAR',tubename)
5198 open (itube,file=tubename,status='old',action='read')
5199 call getenv_loc('IONPAR',ionname)
5200 open (iion,file=ionname,status='old',action='read')
5203 call getenv_loc('ROTPARPDB',rotname_pdb)
5204 open (irotam_pdb,file=rotname_pdb,status='old',action='read')
5207 call getenv_loc('SCPPAR_NUCL',scpname_nucl)
5208 #if defined(WINIFL) || defined(WINPGI)
5209 open (iscpp_nucl,file=scpname_nucl,status='old',readonly,shared)
5210 #elif (defined CRAY) || (defined AIX)
5211 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5213 open (iscpp_nucl,file=scpname_nucl,status='old')
5215 open (iscpp_nucl,file=scpname_nucl,status='old',action='read')
5220 ! 8/9/01 In the newest version SCp interaction constants are read from a file
5221 ! Use -DOLDSCP to use hard-coded constants instead.
5223 call getenv_loc('SCPPAR',scpname)
5224 #if defined(WINIFL) || defined(WINPGI)
5225 open (iscpp,file=scpname,status='old',readonly,shared)
5226 #elif (defined CRAY) || (defined AIX)
5227 open (iscpp,file=scpname,status='old',action='read')
5229 open (iscpp,file=scpname,status='old')
5231 open (iscpp,file=scpname,status='old',action='read')
5234 call getenv_loc('PATTERN',patname)
5235 #if defined(WINIFL) || defined(WINPGI)
5236 open (icbase,file=patname,status='old',readonly,shared)
5237 #elif (defined CRAY) || (defined AIX)
5238 open (icbase,file=patname,status='old',action='read')
5240 open (icbase,file=patname,status='old')
5242 open (icbase,file=patname,status='old',action='read')
5245 ! Open output file only for CG processes
5246 ! print *,"Processor",myrank," fg_rank",fg_rank
5247 if (fg_rank.eq.0) then
5249 if (nodes.eq.1) then
5252 npos = dlog10(dfloat(nodes-1))+1
5254 if (npos.lt.3) npos=3
5255 write (liczba,'(i1)') npos
5256 form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba)) &
5258 write (liczba,form) me
5259 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)// &
5260 liczba(:ilen(liczba))
5261 intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5263 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba)) &
5265 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)// &
5266 liczba(:ilen(liczba))//'.mol2'
5267 statname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5268 liczba(:ilen(liczba))//'.stat'
5270 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot) &
5271 //liczba(:ilen(liczba))//'.stat')
5272 rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba)) &
5275 qname=prefix(:lenpre)//'_'//pot(:lenpot)// &
5276 liczba(:ilen(liczba))//'.const'
5281 outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
5282 intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
5283 pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
5284 mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
5285 statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
5287 call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)// &
5289 rest2name=prefix(:ilen(prefix))//'.rst'
5291 qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
5294 #if defined(AIX) || defined(PGI)
5295 if (me.eq.king .or. .not. out1file) &
5296 open(iout,file=outname,status='unknown')
5298 if (fg_rank.gt.0) then
5299 write (liczba,'(i3.3)') myrank/nfgtasks
5300 write (ll,'(bz,i3.3)') fg_rank
5301 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5306 open(igeom,file=intname,status='unknown',position='append')
5307 open(ipdb,file=pdbname,status='unknown')
5308 open(imol2,file=mol2name,status='unknown')
5309 open(istat,file=statname,status='unknown',position='append')
5311 !1out open(iout,file=outname,status='unknown')
5314 if (me.eq.king .or. .not.out1file) &
5315 open(iout,file=outname,status='unknown')
5317 if (fg_rank.gt.0) then
5318 write (liczba,'(i3.3)') myrank/nfgtasks
5319 write (ll,'(bz,i3.3)') fg_rank
5320 open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,&
5325 open(igeom,file=intname,status='unknown',access='append')
5326 open(ipdb,file=pdbname,status='unknown')
5327 open(imol2,file=mol2name,status='unknown')
5328 open(istat,file=statname,status='unknown',access='append')
5330 !1out open(iout,file=outname,status='unknown')
5333 csa_rbank=prefix(:lenpre)//'.CSA.rbank'
5334 csa_seed=prefix(:lenpre)//'.CSA.seed'
5335 csa_history=prefix(:lenpre)//'.CSA.history'
5336 csa_bank=prefix(:lenpre)//'.CSA.bank'
5337 csa_bank1=prefix(:lenpre)//'.CSA.bank1'
5338 csa_alpha=prefix(:lenpre)//'.CSA.alpha'
5339 csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
5340 !!bankt csa_bankt=prefix(:lenpre)//'.CSA.bankt'
5341 csa_int=prefix(:lenpre)//'.int'
5342 csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
5343 csa_native_int=prefix(:lenpre)//'.CSA.native.int'
5344 csa_in=prefix(:lenpre)//'.CSA.in'
5345 ! print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
5348 write (iout,'(80(1h-))')
5349 write (iout,'(30x,a)') "FILE ASSIGNMENT"
5350 write (iout,'(80(1h-))')
5351 write (iout,*) "Input file : ",&
5352 pref_orig(:ilen(pref_orig))//'.inp'
5353 write (iout,*) "Output file : ",&
5354 outname(:ilen(outname))
5356 write (iout,*) "Sidechain potential file : ",&
5357 sidename(:ilen(sidename))
5359 write (iout,*) "SCp potential file : ",&
5360 scpname(:ilen(scpname))
5362 write (iout,*) "Electrostatic potential file : ",&
5363 elename(:ilen(elename))
5364 write (iout,*) "Cumulant coefficient file : ",&
5365 fouriername(:ilen(fouriername))
5366 write (iout,*) "Torsional parameter file : ",&
5367 torname(:ilen(torname))
5368 write (iout,*) "Double torsional parameter file : ",&
5369 tordname(:ilen(tordname))
5370 write (iout,*) "SCCOR parameter file : ",&
5371 sccorname(:ilen(sccorname))
5372 write (iout,*) "Bond & inertia constant file : ",&
5373 bondname(:ilen(bondname))
5374 write (iout,*) "Bending parameter file : ",&
5375 thetname(:ilen(thetname))
5376 write (iout,*) "Rotamer parameter file : ",&
5377 rotname(:ilen(rotname))
5380 write (iout,*) "Thetpdb parameter file : ",&
5381 thetname_pdb(:ilen(thetname_pdb))
5384 write (iout,*) "Threading database : ",&
5385 patname(:ilen(patname))
5387 write (iout,*)" DIRTMP : ",&
5389 write (iout,'(80(1h-))')
5392 end subroutine openunits
5393 !-----------------------------------------------------------------------------
5396 use geometry_data, only: nres,dc
5398 ! implicit real*8 (a-h,o-z)
5399 ! include 'DIMENSIONS'
5400 ! include 'COMMON.CHAIN'
5401 ! include 'COMMON.IOUNITS'
5402 ! include 'COMMON.MD'
5405 ! real(kind=8) :: var,ene
5407 open(irest2,file=rest2name,status='unknown')
5408 read(irest2,*) totT,EK,potE,totE,t_bath
5411 ! AL 4/17/17: Now reading d_t(0,:) too
5413 read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
5416 ! AL 4/17/17: Now reading d_c(0,:) too
5418 read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
5421 read (irest2,*) iset
5425 end subroutine readrst
5426 !-----------------------------------------------------------------------------
5427 subroutine read_fragments
5431 use control_data, only:out1file
5434 ! implicit real*8 (a-h,o-z)
5435 ! include 'DIMENSIONS'
5439 ! include 'COMMON.SETUP'
5440 ! include 'COMMON.CHAIN'
5441 ! include 'COMMON.IOUNITS'
5442 ! include 'COMMON.MD'
5443 ! include 'COMMON.CONTROL'
5446 ! real(kind=8) :: var,ene
5448 read(inp,*) nset,nfrag,npair,nfrag_back
5450 !el from module energy
5451 ! if(.not.allocated(mset)) allocate(mset(nset)) !(maxprocs/20)
5452 if(.not.allocated(wfrag_back)) then
5453 allocate(wfrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5454 allocate(ifrag_back(3,nfrag_back,nset)) !(3,maxfrag_back,maxprocs/20)
5456 allocate(qinfrag(nfrag,nset),wfrag(nfrag,nset)) !(50,maxprocs/20)
5457 allocate(ifrag(2,nfrag,nset)) !(2,50,maxprocs/20)
5459 allocate(qinpair(npair,nset),wpair(npair,nset)) !(100,maxprocs/20)
5460 allocate(ipair(2,npair,nset)) !(2,100,maxprocs/20)
5463 if(me.eq.king.or..not.out1file) &
5464 write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,&
5465 " nfrag_back",nfrag_back
5467 read(inp,*) mset(iset)
5469 read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),&
5471 if(me.eq.king.or..not.out1file) &
5472 write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),&
5473 ifrag(2,i,iset), qinfrag(i,iset)
5476 read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),&
5478 if(me.eq.king.or..not.out1file) &
5479 write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),&
5480 ipair(2,i,iset), qinpair(i,iset)
5483 read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5484 wfrag_back(3,i,iset),&
5485 ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5486 if(me.eq.king.or..not.out1file) &
5487 write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),&
5488 wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
5492 end subroutine read_fragments
5493 !-----------------------------------------------------------------------------
5495 !-----------------------------------------------------------------------------
5499 ! implicit real*8 (a-h,o-z)
5500 ! include 'DIMENSIONS'
5501 ! include 'COMMON.CSA'
5502 ! include 'COMMON.BANK'
5503 ! include 'COMMON.IOUNITS'
5505 ! integer :: ntf,ik,iw_pdb
5506 ! real(kind=8) :: var,ene
5508 open(icsa_in,file=csa_in,status="old",err=100)
5509 read(icsa_in,*) nconf
5510 read(icsa_in,*) jstart,jend
5511 read(icsa_in,*) nstmax
5512 read(icsa_in,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5513 read(icsa_in,*) nran0,nran1,irr
5514 read(icsa_in,*) nseed
5515 read(icsa_in,*) ntotal,cut1,cut2
5516 read(icsa_in,*) estop
5517 read(icsa_in,*) icmax,irestart
5518 read(icsa_in,*) ntbankm,dele,difcut
5519 read(icsa_in,*) iref,rmscut,pnccut
5520 read(icsa_in,*) ndiff
5527 end subroutine csa_read
5528 !-----------------------------------------------------------------------------
5529 subroutine initial_write
5532 ! implicit real*8 (a-h,o-z)
5533 ! include 'DIMENSIONS'
5534 ! include 'COMMON.CSA'
5535 ! include 'COMMON.BANK'
5536 ! include 'COMMON.IOUNITS'
5538 ! integer :: ntf,ik,iw_pdb
5539 ! real(kind=8) :: var,ene
5541 open(icsa_seed,file=csa_seed,status="unknown")
5542 write(icsa_seed,*) "seed"
5544 #if defined(AIX) || defined(PGI)
5545 open(icsa_history,file=csa_history,status="unknown",&
5548 open(icsa_history,file=csa_history,status="unknown",&
5551 write(icsa_history,*) nconf
5552 write(icsa_history,*) jstart,jend
5553 write(icsa_history,*) nstmax
5554 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5555 write(icsa_history,*) nran0,nran1,irr
5556 write(icsa_history,*) nseed
5557 write(icsa_history,*) ntotal,cut1,cut2
5558 write(icsa_history,*) estop
5559 write(icsa_history,*) icmax,irestart
5560 write(icsa_history,*) ntbankm,dele,difcut
5561 write(icsa_history,*) iref,rmscut,pnccut
5562 write(icsa_history,*) ndiff
5564 write(icsa_history,*)
5567 open(icsa_bank1,file=csa_bank1,status="unknown")
5568 write(icsa_bank1,*) 0
5572 end subroutine initial_write
5573 !-----------------------------------------------------------------------------
5574 subroutine restart_write
5577 ! implicit real*8 (a-h,o-z)
5578 ! include 'DIMENSIONS'
5579 ! include 'COMMON.IOUNITS'
5580 ! include 'COMMON.CSA'
5581 ! include 'COMMON.BANK'
5583 ! integer :: ntf,ik,iw_pdb
5584 ! real(kind=8) :: var,ene
5586 #if defined(AIX) || defined(PGI)
5587 open(icsa_history,file=csa_history,position="append")
5589 open(icsa_history,file=csa_history,access="append")
5591 write(icsa_history,*)
5592 write(icsa_history,*) "This is restart"
5593 write(icsa_history,*)
5594 write(icsa_history,*) nconf
5595 write(icsa_history,*) jstart,jend
5596 write(icsa_history,*) nstmax
5597 write(icsa_history,*) n1,n2,n3,n4,n5,n6,n7,n8,is1,is2
5598 write(icsa_history,*) nran0,nran1,irr
5599 write(icsa_history,*) nseed
5600 write(icsa_history,*) ntotal,cut1,cut2
5601 write(icsa_history,*) estop
5602 write(icsa_history,*) icmax,irestart
5603 write(icsa_history,*) ntbankm,dele,difcut
5604 write(icsa_history,*) iref,rmscut,pnccut
5605 write(icsa_history,*) ndiff
5606 write(icsa_history,*)
5607 write(icsa_history,*) "irestart is: ", irestart
5609 write(icsa_history,*)
5613 end subroutine restart_write
5614 !-----------------------------------------------------------------------------
5616 !-----------------------------------------------------------------------------
5617 subroutine write_pdb(npdb,titelloc,ee)
5619 ! implicit real*8 (a-h,o-z)
5620 ! include 'DIMENSIONS'
5621 ! include 'COMMON.IOUNITS'
5622 character(len=50) :: titelloc1
5623 character*(*) :: titelloc
5624 character(len=3) :: zahl
5625 character(len=5) :: liczba5
5627 integer :: npdb !,ilen
5631 ! real(kind=8) :: var,ene
5635 if (npdb.lt.1000) then
5636 call numstr(npdb,zahl)
5637 open(ipdb,file=prefix(:lenpre)//'@@'//zahl//'.pdb')
5639 if (npdb.lt.10000) then
5640 write(liczba5,'(i1,i4)') 0,npdb
5642 write(liczba5,'(i5)') npdb
5644 open(ipdb,file=prefix(:lenpre)//'@@'//liczba5//'.pdb')
5646 call pdbout(ee,titelloc1,ipdb)
5649 end subroutine write_pdb
5650 !-----------------------------------------------------------------------------
5652 !-----------------------------------------------------------------------------
5653 subroutine write_thread_summary
5654 ! Thread the sequence through a database of known structures
5655 use control_data, only: refstr
5657 use energy_data, only: n_ene_comp
5659 ! implicit real*8 (a-h,o-z)
5660 ! include 'DIMENSIONS'
5662 use MPI_data !include 'COMMON.INFO'
5665 ! include 'COMMON.CONTROL'
5666 ! include 'COMMON.CHAIN'
5667 ! include 'COMMON.DBASE'
5668 ! include 'COMMON.INTERACT'
5669 ! include 'COMMON.VAR'
5670 ! include 'COMMON.THREAD'
5671 ! include 'COMMON.FFIELD'
5672 ! include 'COMMON.SBRIDGE'
5673 ! include 'COMMON.HEADER'
5674 ! include 'COMMON.NAMES'
5675 ! include 'COMMON.IOUNITS'
5676 ! include 'COMMON.TIME1'
5678 integer,dimension(maxthread) :: ip
5679 real(kind=8),dimension(0:n_ene) :: energia
5681 integer :: i,j,ii,jj,ipj,ik,kk,ist
5682 real(kind=8) :: enet,etot,rmsnat,rms,frac,frac_nn
5684 write (iout,'(30x,a/)') &
5685 ' *********** Summary threading statistics ************'
5686 write (iout,'(a)') 'Initial energies:'
5687 write (iout,'(a4,2x,a12,14a14,3a8)') &
5688 'No','seq',(ename(print_order(i)),i=1,nprint_ene),'ETOT',&
5689 'RMSnat','NatCONT','NNCONT','RMS'
5690 ! Energy sort patterns
5695 enet=ener(n_ene-1,ip(i))
5698 if (ener(n_ene-1,ip(j)).lt.enet) then
5700 enet=ener(n_ene-1,ip(j))
5712 ist=nres_base(2,ii)+ipatt(2,i)
5714 energia(i)=ener0(kk,i)
5716 etot=ener0(n_ene_comp+1,i)
5717 rmsnat=ener0(n_ene_comp+2,i)
5718 rms=ener0(n_ene_comp+3,i)
5719 frac=ener0(n_ene_comp+4,i)
5720 frac_nn=ener0(n_ene_comp+5,i)
5723 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5724 i,str_nam(ii),ist+1,&
5725 (energia(print_order(kk)),kk=1,nprint_ene),&
5726 etot,rmsnat,frac,frac_nn,rms
5728 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3)') &
5729 i,str_nam(ii),ist+1,&
5730 (energia(print_order(kk)),kk=1,nprint_ene),etot
5733 write (iout,'(//a)') 'Final energies:'
5734 write (iout,'(a4,2x,a12,17a14,3a8)') &
5735 'No','seq',(ename(print_order(kk)),kk=1,nprint_ene),'ETOT',&
5736 'RMSnat','NatCONT','NNCONT','RMS'
5740 ist=nres_base(2,ii)+ipatt(2,i)
5742 energia(kk)=ener(kk,ik)
5744 etot=ener(n_ene_comp+1,i)
5745 rmsnat=ener(n_ene_comp+2,i)
5746 rms=ener(n_ene_comp+3,i)
5747 frac=ener(n_ene_comp+4,i)
5748 frac_nn=ener(n_ene_comp+5,i)
5749 write (iout,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5750 i,str_nam(ii),ist+1,&
5751 (energia(print_order(kk)),kk=1,nprint_ene),&
5752 etot,rmsnat,frac,frac_nn,rms
5754 write (iout,'(/a/)') 'IEXAM array:'
5755 write (iout,'(i5)') nexcl
5757 write (iout,'(2i5)') iexam(1,i),iexam(2,i)
5759 write (iout,'(/a,1pe14.4/a,1pe14.4/)') &
5760 'Max. time for threading step ',max_time_for_thread,&
5761 'Average time for threading step: ',ave_time_for_thread
5763 end subroutine write_thread_summary
5764 !-----------------------------------------------------------------------------
5765 subroutine write_stat_thread(ithread,ipattern,ist)
5767 use energy_data, only: n_ene_comp
5769 ! implicit real*8 (a-h,o-z)
5770 ! include "DIMENSIONS"
5771 ! include "COMMON.CONTROL"
5772 ! include "COMMON.IOUNITS"
5773 ! include "COMMON.THREAD"
5774 ! include "COMMON.FFIELD"
5775 ! include "COMMON.DBASE"
5776 ! include "COMMON.NAMES"
5777 real(kind=8),dimension(0:n_ene) :: energia
5779 integer :: ithread,ipattern,ist,i
5780 real(kind=8) :: etot,rmsnat,rms,frac,frac_nn
5782 #if defined(AIX) || defined(PGI)
5783 open(istat,file=statname,position='append')
5785 open(istat,file=statname,access='append')
5788 energia(i)=ener(i,ithread)
5790 etot=ener(n_ene_comp+1,ithread)
5791 rmsnat=ener(n_ene_comp+2,ithread)
5792 rms=ener(n_ene_comp+3,ithread)
5793 frac=ener(n_ene_comp+4,ithread)
5794 frac_nn=ener(n_ene_comp+5,ithread)
5795 write (istat,'(i4,2x,a8,i4,14(1pe14.5),0pf8.3,f8.5,f8.5,f8.3)') &
5796 ithread,str_nam(ipattern),ist+1,&
5797 (energia(print_order(i)),i=1,nprint_ene),&
5798 etot,rmsnat,frac,frac_nn,rms
5801 end subroutine write_stat_thread
5802 !-----------------------------------------------------------------------------
5804 !-----------------------------------------------------------------------------
5805 end module io_config